{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test scipy.signal.correlate on some atl06 data from foundation ice stream"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy, sys, os, pyproj, glob, re, h5py\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy.signal import correlate\n",
    "from astropy.time import Time\n",
    "\n",
    "%matplotlib widget\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test scipy.signal.correlate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate some test data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "20608b8b8ba8467ab8eb991681d789cc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'index')"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dx = 0.1\n",
    "x = np.arange(0,10,dx)\n",
    "y = np.zeros(np.shape(x))\n",
    "ix0 = 30\n",
    "ix1 = 30 + 15\n",
    "y[ix0:ix1] = 1\n",
    "\n",
    "fig,axs = plt.subplots(1,2)\n",
    "axs[0].plot(x,y,'k')\n",
    "axs[0].set_xlabel('distance (m)')\n",
    "axs[0].set_ylabel('value')\n",
    "axs[1].plot(np.arange(len(x)), y,'k')\n",
    "axs[1].set_xlabel('index')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make a signal to correlate with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "40886133dcaf49e9b90be0f113b75597",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0.98, 'black = original, blue = shifted')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "imposed_offset = int(14/dx) # 14 meters, in units of samples\n",
    "\n",
    "x_noise = np.arange(0,50,dx) # make the vector we're comparing with much longer\n",
    "y_noise = np.zeros(np.shape(x_noise))\n",
    "y_noise[ix0 + imposed_offset : ix1 + imposed_offset] = 1\n",
    "\n",
    "# uncomment the line below to add noise\n",
    "# y_noise = y_noise * np.random.random(np.shape(y_noise))\n",
    "\n",
    "fig,axs = plt.subplots(1,2)\n",
    "\n",
    "axs[0].plot(x,y,'k')\n",
    "axs[0].set_xlabel('distance (m)')\n",
    "axs[0].set_ylabel('value')\n",
    "axs[1].plot(np.arange(len(x)), y, 'k')\n",
    "axs[1].set_xlabel('index')\n",
    "\n",
    "axs[0].plot(x_noise,y_noise, 'b')\n",
    "axs[0].set_xlabel('distance (m)')\n",
    "axs[0].set_ylabel('value')\n",
    "axs[1].plot(np.arange(len(x_noise)), y_noise,'b')\n",
    "axs[1].set_xlabel('index')\n",
    "\n",
    "fig.suptitle('black = original, blue = shifted')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try scipy.signal.correlate:\n",
    "\n",
    "mode ='full' returns the entire cross correlation; could be 'valid' to return only non- zero-padded part\n",
    "\n",
    "method = direct (not fft)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "corr = correlate(y_noise,y, mode = 'full', method = 'direct') \n",
    "norm_val = np.sqrt(np.sum(y_noise**2)*np.sum(y**2))\n",
    "corr = corr / norm_val"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What are the dimensions of corr?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "corr:  (599,)\n",
      "x:  (100,)\n",
      "x:  (500,)\n"
     ]
    }
   ],
   "source": [
    "print('corr: ', np.shape(corr))\n",
    "print('x: ', np.shape(x))\n",
    "print('x: ', np.shape(x_noise))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2cb19916a8434cc5915af29ca701d26b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0.98, 'Shift  140  samples, or  14.0  m to line up signals')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# lagvec = np.arange(0,len(x_noise) - len(x) + 1)\n",
    "lagvec = np.arange( -(len(x) - 1), len(x_noise), 1)\n",
    "shift_vec = lagvec * dx\n",
    "\n",
    "ix_peak = np.arange(len(corr))[corr == np.nanmax(corr)][0]\n",
    "best_lag = lagvec[ix_peak]\n",
    "best_shift = shift_vec[ix_peak]\n",
    "\n",
    "fig,axs = plt.subplots(3,1)\n",
    "\n",
    "axs[0].plot(lagvec,corr)\n",
    "axs[0].plot(lagvec[ix_peak],corr[ix_peak], 'r*')\n",
    "axs[0].set_xlabel('lag (samples)')\n",
    "axs[0].set_ylabel('correlation coefficient')\n",
    "\n",
    "axs[1].plot(shift_vec,corr)\n",
    "axs[1].plot(shift_vec[ix_peak],corr[ix_peak], 'r*')\n",
    "axs[1].set_xlabel('shift (m)')\n",
    "axs[1].set_ylabel('correlation coefficient')\n",
    "\n",
    "axs[2].plot(x + best_shift, y,'k')\n",
    "axs[2].plot(x_noise, y_noise, 'b--')\n",
    "axs[2].set_xlabel('shift (m)')\n",
    "\n",
    "\n",
    "fig.suptitle(' '.join(['Shift ', str(best_lag), ' samples, or ', str(best_shift), ' m to line up signals']))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let's try with some ATL06 data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load some repeat data:\n",
    "\n",
    "\n",
    "import readers, etc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#! cd ..; [ -d pointCollection ] || git clone https://www.github.com/smithB/pointCollection.git\n",
    "#sys.path.append(os.path.join(os.getcwd(), '..'))\n",
    "\n",
    "#!python3 -m pip install --user git+https://github.com/tsutterley/pointCollection.git@pip\n",
    "import pointCollection as pc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "moa_datapath = '/srv/tutorial-data/land_ice_applications/'\n",
    "datapath = '/home/jovyan/shared/surface_velocity/FIS_ATL06/'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# example hf5 file, if you need to look at the fields\n",
    "# datapath='/home/jovyan/shared/surface_velocity/FIS_ATL06_small/processed_ATL06_20191129105346_09700511_003_01.h5'\n",
    "# !h5ls -r /home/jovyan/shared/surface_velocity/FIS_ATL06_small/processed_ATL06_20191129105346_09700511_003_01.h5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Geographic setting : Foundation Ice Stream"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/jovyan/.local/lib/python3.7/site-packages/pointCollection/__init__.py\n"
     ]
    }
   ],
   "source": [
    "print(pc.__file__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-76.  -74.5 -74.5 -76.  -76. ]\n",
      "[ -98.  -98. -102. -102.  -98.]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "45c67f55989d4ce0a6571b5de0570bca",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'cmap': 'gray', 'clim': [14000, 17000], 'extent': array([-1676950., -1495950.,  -352175.,  -213175.]), 'origin': 'lower'}\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Mosaic of Antarctica for Pine Island Glacier')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# something wrong with pointCollection\n",
    "\n",
    "spatial_extent = np.array([-102, -76, -98, -74.5])\n",
    "lat=spatial_extent[[1, 3, 3, 1, 1]]\n",
    "lon=spatial_extent[[2, 2, 0, 0, 2]]\n",
    "print(lat)\n",
    "print(lon)\n",
    "# project the coordinates to Antarctic polar stereographic\n",
    "xy=np.array(pyproj.Proj(3031)(lon, lat))\n",
    "# get the bounds of the projected coordinates \n",
    "XR=[np.nanmin(xy[0,:]), np.nanmax(xy[0,:])]\n",
    "YR=[np.nanmin(xy[1,:]), np.nanmax(xy[1,:])]\n",
    "MOA=pc.grid.data().from_geotif(os.path.join(moa_datapath, 'MOA','moa_2009_1km.tif'), bounds=[XR, YR])\n",
    "\n",
    "# show the mosaic:\n",
    "plt.figure()\n",
    "MOA.show(cmap='gray', clim=[14000, 17000])\n",
    "plt.plot(xy[0,:], xy[1,:])\n",
    "plt.title('Mosaic of Antarctica for Pine Island Glacier')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load repeat track data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ATL06 reader"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def atl06_to_dict(filename, beam, field_dict=None, index=None, epsg=None):\n",
    "    \"\"\"\n",
    "        Read selected datasets from an ATL06 file\n",
    "\n",
    "        Input arguments:\n",
    "            filename: ATl06 file to read\n",
    "            beam: a string specifying which beam is to be read (ex: gt1l, gt1r, gt2l, etc)\n",
    "            field_dict: A dictinary describing the fields to be read\n",
    "                    keys give the group names to be read, \n",
    "                    entries are lists of datasets within the groups\n",
    "            index: which entries in each field to read\n",
    "            epsg: an EPSG code specifying a projection (see www.epsg.org).  Good choices are:\n",
    "                for Greenland, 3413 (polar stereographic projection, with Greenland along the Y axis)\n",
    "                for Antarctica, 3031 (polar stereographic projection, centered on the Pouth Pole)\n",
    "        Output argument:\n",
    "            D6: dictionary containing ATL06 data.  Each dataset in \n",
    "                dataset_dict has its own entry in D6.  Each dataset \n",
    "                in D6 contains a numpy array containing the \n",
    "                data\n",
    "    \"\"\"\n",
    "    if field_dict is None:\n",
    "        field_dict={None:['latitude','longitude','h_li', 'atl06_quality_summary'],\\\n",
    "                    'ground_track':['x_atc','y_atc'],\\\n",
    "                    'fit_statistics':['dh_fit_dx', 'dh_fit_dy']}\n",
    "    D={}\n",
    "    # below: file_re = regular expression, it will pull apart the regular expression to get the information from the filename\n",
    "    file_re=re.compile('ATL06_(?P<date>\\d+)_(?P<rgt>\\d\\d\\d\\d)(?P<cycle>\\d\\d)(?P<region>\\d\\d)_(?P<release>\\d\\d\\d)_(?P<version>\\d\\d).h5')\n",
    "    with h5py.File(filename,'r') as h5f:\n",
    "        for key in field_dict:\n",
    "            for ds in field_dict[key]:\n",
    "                if key is not None:\n",
    "                    ds_name=beam+'/land_ice_segments/'+key+'/'+ds\n",
    "                else:\n",
    "                    ds_name=beam+'/land_ice_segments/'+ds\n",
    "                if index is not None:\n",
    "                    D[ds]=np.array(h5f[ds_name][index])\n",
    "                else:\n",
    "                    D[ds]=np.array(h5f[ds_name])\n",
    "                if '_FillValue' in h5f[ds_name].attrs:\n",
    "                    bad_vals=D[ds]==h5f[ds_name].attrs['_FillValue']\n",
    "                    D[ds]=D[ds].astype(float)\n",
    "                    D[ds][bad_vals]=np.NaN\n",
    "        D['data_start_utc'] = h5f['/ancillary_data/data_start_utc'][:]\n",
    "        D['delta_time'] = h5f['/' + beam + '/land_ice_segments/delta_time'][:]\n",
    "        D['segment_id'] = h5f['/' + beam + '/land_ice_segments/segment_id'][:]\n",
    "    if epsg is not None:\n",
    "        xy=np.array(pyproj.proj.Proj(epsg)(D['longitude'], D['latitude']))\n",
    "        D['x']=xy[0,:].reshape(D['latitude'].shape)\n",
    "        D['y']=xy[1,:].reshape(D['latitude'].shape)\n",
    "    temp=file_re.search(filename)\n",
    "    D['rgt']=int(temp['rgt'])\n",
    "    D['cycle']=int(temp['cycle'])\n",
    "    D['beam']=beam\n",
    "    return D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read in files; this next cell took ~1 minute early in the morning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "file /home/jovyan/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190430122344_04920311_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /home/jovyan/shared/surface_velocity/FIS_ATL06/processed_ATL06_20181030210407_04920111_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /home/jovyan/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190730080323_04920411_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /home/jovyan/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190220230230_08320211_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /home/jovyan/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190312235510_11380211_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /home/jovyan/shared/surface_velocity/FIS_ATL06/processed_ATL06_20181108184743_06280111_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /home/jovyan/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190228224553_09540211_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /home/jovyan/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190623094402_13150311_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /home/jovyan/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190506112405_05830311_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /home/jovyan/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190611193446_11380311_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /home/jovyan/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190620171809_12740311_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /home/jovyan/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190823150456_08630411_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /home/jovyan/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190529105941_09340311_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /home/jovyan/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190308130329_10700211_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /home/jovyan/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190824143917_08780411_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /home/jovyan/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190812071246_06900411_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /home/jovyan/shared/surface_velocity/FIS_ATL06/processed_ATL06_20181118033101_07710111_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /home/jovyan/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190822060451_08420411_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /home/jovyan/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190607194306_10770311_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "file /home/jovyan/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190129164404_04920211_003_01.h5 encountered error 'Unable to open object (component not found)'\n",
      "read 613 data files of which 20 gave errors\n"
     ]
    }
   ],
   "source": [
    "# find all the files in the directory:\n",
    "# ATL06_files=glob.glob(os.path.join(datapath, 'PIG_ATL06', '*.h5'))\n",
    "\n",
    "ATL06_files=glob.glob(os.path.join(datapath, '*.h5'))\n",
    "\n",
    "D_dict={}\n",
    "error_count=0\n",
    "for file in ATL06_files:\n",
    "    try:\n",
    "        D_dict[file]=atl06_to_dict(file, '/gt2l', index=slice(0, -1, 25), epsg=3031)\n",
    "    except KeyError as e:\n",
    "        print(f'file {file} encountered error {e}')\n",
    "        error_count += 1\n",
    "print(f\"read {len(D_dict)} data files of which {error_count} gave errors\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot ground tracks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1db02f4484bc4dc09d9441bffa8a6015",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'cmap': 'gray', 'clim': [14000, 17000], 'extent': array([-1676950., -1495950.,  -352175.,  -213175.]), 'origin': 'lower'}\n",
      "{'cmap': 'gray', 'clim': [14000, 17000], 'extent': array([-1676950., -1495950.,  -352175.,  -213175.]), 'origin': 'lower'}\n"
     ]
    }
   ],
   "source": [
    "plt.figure(figsize=[8,8])\n",
    "hax0=plt.gcf().add_subplot(211, aspect='equal')\n",
    "MOA.show(ax=hax0, cmap='gray', clim=[14000, 17000]);\n",
    "hax1=plt.gcf().add_subplot(212, aspect='equal', sharex=hax0, sharey=hax0)\n",
    "MOA.show(ax=hax1, cmap='gray', clim=[14000, 17000]);\n",
    "for fname, Di in D_dict.items():\n",
    "    cycle=Di['cycle']\n",
    "    if cycle <= 2:\n",
    "        ax=hax0\n",
    "    else:\n",
    "        ax=hax1\n",
    "    #print(fname)\n",
    "    #print(f'\\t{rgt}, {cycle}, {region}')\n",
    "    ax.plot(Di['x'], Di['y'])\n",
    "    if True:\n",
    "        try:\n",
    "            if cycle  < 3:\n",
    "                ax.text(Di['x'][0], Di['y'][0], f\"rgt={Di['rgt']}, cyc={cycle}\", clip_on=True)\n",
    "            elif cycle==3:\n",
    "                ax.text(Di['x'][0], Di['y'][0], f\"rgt={Di['rgt']}, cyc={cycle}+\", clip_on=True)\n",
    "        except IndexError:\n",
    "            pass\n",
    "hax0.set_title('cycles 1 and 2');\n",
    "hax1.set_title('cycle 3+');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Map view elevations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b0e54914932d499f870f35862acc97c8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "map_fig=plt.figure()\n",
    "map_ax=map_fig.add_subplot(111)\n",
    "# MOA.show(ax=map_ax, cmap='gray', clim=[14000, 17000])\n",
    "for fname, Di in D_dict.items():\n",
    "    # select elevations with good quality_summary\n",
    "    good=Di['atl06_quality_summary']==0\n",
    "    ms=map_ax.scatter( Di['x'][good], Di['y'][good],  2, c=Di['h_li'][good], \\\n",
    "                  vmin=0, vmax=1000, label=fname)\n",
    "map_ax._aspect='equal'\n",
    "plt.colorbar(ms, label='elevation');\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Repeat track elevation profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ben Smiths's code to plot the individual segments:\n",
    "def plot_segs(D6, ind=None, **kwargs):\n",
    "    \"\"\"\n",
    "    Plot a sloping line for each ATL06 segment\n",
    "    \"\"\"\n",
    "    if ind is None:\n",
    "        ind=np.ones_like(D6['h_li'], dtype=bool)\n",
    "    #define the heights of the segment endpoints.  Leave a row of NaNs so that the endpoints don't get joined\n",
    "    h_ep=np.zeros([3, D6['h_li'][ind].size])+np.NaN\n",
    "    h_ep[0, :]=D6['h_li'][ind]-D6['dh_fit_dx'][ind]*20\n",
    "    h_ep[1, :]=D6['h_li'][ind]+D6['dh_fit_dx'][ind]*20\n",
    "    # define the x coordinates of the segment endpoints\n",
    "    x_ep=np.zeros([3,D6['h_li'][ind].size])+np.NaN\n",
    "    x_ep[0, :]=D6['x_atc'][ind]-20\n",
    "    x_ep[1, :]=D6['x_atc'][ind]+20\n",
    "\n",
    "    plt.plot(x_ep.T.ravel(), h_ep.T.ravel(), **kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A revised code to plot the elevations of segment midpoints (h_li):\n",
    "def plot_elevation(D6, ind=None, **kwargs):\n",
    "    \"\"\"\n",
    "    Plot midpoint elevation for each ATL06 segment\n",
    "    \"\"\"\n",
    "    if ind is None:\n",
    "        ind=np.ones_like(D6['h_li'], dtype=bool)\n",
    "    # pull out heights of segment midpoints\n",
    "    h_li = D6['h_li'][ind]\n",
    "    # pull out along track x coordinates of segment midpoints\n",
    "    x_atc = D6['x_atc'][ind]\n",
    "\n",
    "    plt.plot(x_atc, h_li, **kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "492aed7c78fb422d8daaa41f1e2f04ec",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "D_2l={}\n",
    "D_2r={}\n",
    "\n",
    "# specify the rgt here:\n",
    "rgt=\"0027\"\n",
    "rgt=\"0848\" #Ben's suggestion\n",
    "\n",
    "# iterate over the repeat cycles\n",
    "for cycle in ['03','04','05','06','07']:\n",
    "    for filename in glob.glob(os.path.join(datapath, f'*ATL06_*_{rgt}{cycle}*_003*.h5')):\n",
    "        try:\n",
    "            # read the left-beam data\n",
    "            D_2l[filename]=atl06_to_dict(filename,'/gt2l', index=None, epsg=3031)\n",
    "            # read the right-beam data\n",
    "            D_2r[filename]=atl06_to_dict(filename,'/gt2r', index=None, epsg=3031)\n",
    "            # plot the locations in the previous plot\n",
    "            map_ax.plot(D_2r[filename]['x'], D_2r[filename]['y'],'k');  \n",
    "            map_ax.plot(D_2l[filename]['x'], D_2l[filename]['y'],'k');\n",
    "        except Exception as e:\n",
    "            print(f'filename={filename}, exception={e}')\n",
    "\n",
    "plt.figure();\n",
    "for filename, Di in D_2l.items():\n",
    "    #Plot only points that have ATL06_quality_summary==0 (good points)\n",
    "    hl=plot_elevation(Di, ind=Di['atl06_quality_summary']==0, label=f\"cycle={Di['cycle']}\")\n",
    "    #hl=plt.plot(Di['x_atc'][Di['atl06_quality_summary']==0], Di['h_li'][Di['atl06_quality_summary']==0], '.', label=f\"cycle={Di['cycle']}\")\n",
    "    \n",
    "plt.legend()\n",
    "plt.xlabel('x_atc')\n",
    "plt.ylabel('elevation');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pull out a segment and cross correlate: \n",
    "\n",
    "Let's try x_atc = 2.935e7 thru 2.93e7 (just from looking through data)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deZwcVbnw8d9T3T3TM5nJZJKZLGSbQAIhAULiyBaEIC6AKNz3qqBXREUjIiCKV/HeCwa4CoK4ICgGQeFKRBZBVlnDDgkTyErIPtmXmSSzZZZe6nn/6Jq1e/aenqWf7+fTTPWpU3VOnTRPnz5VdUpUFWOMMenB6e8KGGOMSR0L+sYYk0Ys6BtjTBqxoG+MMWnEgr4xxqQRf39XoDMFBQVaVFTU39UwxphBZdmyZeWqWtg2fcAH/aKiIkpKSvq7GsYYkzKb16wjIsqRM6b3eB8isjVRug3vGGPMAOK6Lvc//DcWPfQg4Yb6pO9/wPf0jTEmnTiOw4zIBByEqoMNjBobTOr+LegbY8wAsn/XAXI1i1yFUWPzkr5/C/rGGDOA3PPHO6kNRJkRHtsn++90TF9E7hWRfSKyukXaSBF5QUQ2eH/zW6z7iYhsFJF1IvLpFukfEZFV3rrbRUSSfzjGGDO41UoUgGMDdX2y/66cyP0LcFabtGuAl1R1GvCS9x4RmQFcCMz0tvm9iPi8bf4AzAemea+2+zTGmLRUXdXATb94nSf+7xEARrvDCfrCfVJWp0FfVV8DDrRJPg+4z1u+Dzi/RfqDqtqgqluAjcAJIjIOGK6qb2tsWs/7W2xjjDFpbe3iUk6u3MN7m2IDKqdoGWO+enaflNXTMf0xqrobQFV3i8hoL3088E6LfDu8tLC33DY9IRGZT+xXAZMmTephFY0xZnCocQ/xWsaapvfH33hTn5WV7Ov0E43TawfpCanqQlUtVtXiwsK4G8qMMWZIie56uGl5wYIFfVpWT4P+Xm/IBu/vPi99BzCxRb4JwC4vfUKCdGOMSXuH9sRGNM5p6PkduF3V06D/BHCxt3wx8M8W6ReKSKaITCF2wnapNxRULSIneVftfLXFNsYYMzQ01MCy+6CLTySMRCIsWLCA5b69AAz/0ml9WTuga5ds/g14GzhKRHaIyCXAzcAnRWQD8EnvPaq6BngI+AD4F/BdVY16u/oO8CdiJ3c3Ac8m+ViMMaZ/vbgAnrwSNi/uNGs4HOYHTyzAOa0Uny8EwPTjDuvjCnbhRK6qfqmdVWe2k/9nwM8SpJcAx3SrdsYYM4hs3VHHY/o9zl61nIy1a5nymcugnVuSrvjDHTx+7OcBuG3Cb5hadElK6mh35BpjTCfKdlXj8zuMHD2sw3z7dx9NhZSzdOlqLsp4lA8Kiplx0omt8ix/bj2vrVnH47PPaEo7/RNXMWHi8X1S97Zslk1jjOlE1e3vsv/Xb3eab5jAxfXzgFoAKravi8tTsmYdP5/dfL3L7nmzUhbwwXr6xhjToXA4wsOZb3NEdCzTOsnrC7gE6nxkS2xK5KiT2bTunt99lZsO+wY1LQL+njNSF+wbWdA3xpgOlK9fwYzIBAo1j3A4SiDgazevOi4ARzmx55dURbNhQWymzEuAS/b/k7GnvwrAxpnj+rbi7bDhHWOM6cCIJ69iTvRwJrqj2Lu7ptU6143w9LLLKd0XC+REYydtH4+exe2R8/nMB19plf/s437Ppa8+zvo5ReSMHpOS+rdlQd8YYzqwv/DzTcsBbWi9rqaUYOWzLH3/Mg7U7oCosjiwmkfdkznNWdkq79jTX+WjK9ezYMEChueNSEndE7HhHWOM6YBbEWjqHQfaPL6wYt8K/sAVHBX4kFHvnE7g+HFIydG8H7ymKc9BJ5ejP/YUF7/1NFd/7/sprHliFvSNMaYDzsHmE6+V9T5GtlhXc8jlDZnHG8zj4/oCoZzdfMVZ1Wr7409+hO2nzsR3+nE4Tv8PrvR/DYwxZgCpev9Vtl3zKlt+8QwAft9S9ks1S/wb2L+3rCnf09v34P7xbk7VVzhO3+c/5FGOeat5+oXzi25j3MdeZusnTyIQCAyIgA/W0zfGmCbb//grntnhY3fwILOrixh7sIZ3okfQ4NvJWv9OAlsygOO5/Pl/8UhgLHzhFu6OXsQJ75cz5lBtq33NKt3Bj4pm9M+BdEC0ixMD9Zfi4mItKSnp72oYY4a4hrow9950O3ud6qY0IX4O+Hq/n7/MPZfi8hU8tebKuP3sjZ7GmBuf7NvKdoGILFPV4rbpvfq9ISKl3nNvl4tIiZfW7efnGmNMf1v1l8ebAv43v/hZoDHgKyuPmcSFwcc4Lft1bo7cyp5XT48L+G8ecRjb65+kQi5IbcW7KRnDO2eoanmL943Pz71ZRK7x3v+4zfNzDwNeFJEjW8zCaYwx/eLla18i6kb4KGM447T9ZE09mhH8lbN0MdNlC8SeYsh0Sltt99ywE7jh2G+xOeMILuUOdNPLnBd0U38A3dAXY/rnAfO85fuAV4Af0+L5ucAWEdkInEBs2mZjjOkXrqu85nsdfDBW95L99iJ4+1augrhn/j2ccSXnfGM+w0aPB8fhzIYwW3/2S2qP/hdOuXAAqK/9sB+Oout6G/QVeF5EFPijqi6k+8/PjWPPyDXG9IWa6mp+edttTe+vvfZaJAp+DfM/ckfCB7tW+oZx+il/ZfmZH+cLbdb5MwPMv+EnQPNjDoOzTu6j2idHb4P+XFXd5QX2F0Sko6+4Lj8n1/vyWAixE7m9rKMxJk099cSTlLy3rN31f77hUr4pD/E/baLTYR97iZ0f/wgiQh6wvAtlff28ebzxwhNMnXdpr+rc13oV9FV1l/d3n4g8Rmy4Zq+IjPN6+V15fq4xxiRVeVkZd9x5JwB3nX5+q3Vf3vkkv9r4y9ibNsH+kunX8/SYeWyaOxNp5+En7Zk8ex6TZ8/raZVTpsdBX0SGAY6qVnvLnwJuoPn5uTcT//zcRSLyK2IncqcBS3tRd2NMmlHVToOxui533Hknd51+Prnhav5n/R+4fPeD7eYv843goyc/xI21W7jncxcmu8oDTm96+mOAx7x/AD+wSFX/JSLvAg95z9LdBrFhMFVdIyKNz8+N0Pr5ucYYE+fDD9bz4EOLvHdKy675tddeC4f2sXTRLTy3J68pPeTz8cacWex59fRO93+9XIFG/XztjWe5yBuTH+rs5ixjzICh0Qg1yx4ltPltatc+y0TZk9wCFlS2ertt+w7GjhlNRkZGcssZANq7OcumYTDGpFR9NErRa6u4evOf+c/tf2m1ToBcb3lU94bUE2sT5NuaNHFCEgoZXCzoG2OSSxX30H4OLllE/uvX4bS5SC8IdLX/viMwmt2ZBdw18QJW5c1mW2ZeXJ4R4Sr8hxqYVLuL32+9ia8f87+8cO6/43eS8a0x9FjQN8Z0W30owv63HyKy8lEm738xbr0DjOrG/m7O/BZvjz6eyfv3MrK2mk0F47jgvM9yRH4+xcOC3JPg5O3zZRUcHgwwNXdYi9SLWdzto0kvFvSNMVB3EDJywfHRUHeILUteYPqrX2s3e5B27qzsxKwTH2HOhg1cPmcOHzn1pKb0b+7dQ/0f7mp6f/9353d6lc6nCvvv6VODmZ3INSaNRMNh5GejcUje/DD/NfFynh43j71ZhQD8cHiIH8wuHjDzx6crO5FrTJpRVZ7atom1O9cje9fxnysW4OvG9quCR/BW/vHcNuVrHLt9I4fv2cOwcIgrr/ge+SNHNPXEf+69zOBgQd+YQaaqtpbXtq/mjfJyInXVFFSV8oDM5sQDK/l45Vtsl7H8uPIBBBjjTOcwdfmIrm+1j1smfJmdwakcP7qI08cfRaQhzMQjjyQrp3l8/Fjv9W0ATkvdAZo+ZUHfmC5SVRrq6glmZxFpqOdQTTUNNTUcqqygpqqSujqXvMKxTJo6iZ3bt/HO5qXsaKhimAMjauoZl38EtXKIkh0fUB9VshsayIvUEFYffokSdgK44qMmIxvXJ+BCbvgQ+eFqstwGhkXryI3UclZ0KecCk5lCoVYwmgqukcTDtGH8RFDe4yjmsI6tM69i8heu50epbTozgFjQN0ml4QYi1ZXUaiYN9SGq925Gxc+huhDb9q0nL+LQEPazp2IL4YYK6lSpdYKgPoZFK3DVJerLZHj2cAqyCtkb3Y9U7aI2KgSidYAPkSgRESKOH/ABDiF/BiNDIUIiVDk+ciJ1DIvWokAow089GWSFHQrrq4i6Qr3j51BOgJqMLKJ+h/LMIOoPkOEPEig/wLSqUgr1IAcCuRzyZbMvcyR54WqOritlTmg9eVJL/MWDEFVBEaaKy1QvzVXBaRGUP9eDdq3TDCoZRrVkN6WV+fLZIuNpiGZx0D+KkW49eZFadmXP4uhjZpFfNIkTj5qJ39/8v/nkHpRthpYhG/T/9qfbKa3eTjBSi4OLtLhWWAEXUG9M0lEl9r+qIqoogiuCoPhQkNg6VPERjW0vjre9oirEZpcWRGP7cGJ7iZWttCq/cbmxLEFxVJvq6agC2nTDuWqsPBcBBEdcVGNrfbj4icSuhdbG42tz1YM038CujQmxCoC0mXMqdpixeilkaogMDZOpYYLecoAIAY0SpIFRWoUiRHEYRj1BCRMAhqlDDsroFsFuVnf+ATu+p6Z76hKkuUBV+5uE1ccuRpEXPsRwanFqYsexRcfyQkYx24PjqHaG0eAECDsB1Bdrq8L6g2S6DdT5g+Q11JMdqceHy87gaHCiBCMhpD4LXyhIvQQJ4WOkRMkUP44EEHFRXwP+YBaOL5tAIEh2Xh7BrFyyRuQxetrhhMdPIOD38fEkNpFJH0M26K+gmk9Wl3Cm+15/VyUhV2OBsvEqigi+xq8I7+uiMRQ3BvcoDrEvk8btGoNtpGltcm9GUYR6yaCBQOwlGTRIgFoJEsFHyAmw3ze86Yuq3smgTjKpc4KMjFYDSlkgv+nI/FGXkOMnJ1pLHUHqyCZTXbJpwE+YCicPQfC7UdQJ4/qVjJBSyQgCEqXOycT1uWjUISMEGeridyK4PpBohMqcTHwRl/y6Bup8WYScIC6KP6IEqSc8zGVrQSH+UIScunqGVUXJDtXjJ0JWyCVD6slw6oiSi2SMwYlmURdR/P4AhYePY+qck5hy1BSm+Bwi4Qh76ipxG+oYkT2MrGAOfsff7ZkZjUm1IRv0z544k9V19WwIn4g6flQEcQQfDo44BAQyAhn4HAf8PkCIOg6okh3wk5WRSchVIq5LMDOHzGAeuXn55ObmEfDyqThk+n34/T4iKojj4PP78Dk+HL8Pvz8TJ5CB44No1AEEn1/wiYPjd/D7fbEwLZDZj8FCVVGlqS4WuDrnD/iZEOjO7UfGDAx2nb4xxgxB7V2nP+CDvoiUAVv7ux4DSAFQ3mmu9GJt0pq1R7x0bJPJqlrYNnHAB33TmoiUJPr2TmfWJq1Ze8SzNmlm90kbY0wasaBvjDFpxIL+4LOwvyswAFmbtGbtEc/axGNj+sYYk0asp2+MMWnEgv4AICITRWSxiKwVkTUi8r0EefJF5DERWSkiS0XkmK5uOxj1sk2C3vsV3rbXp/4Ikq83bdJivU9E3heRp1JX877T2zYRkVIRWSUiy0UkPW4Iit2Naa/+fAHjgDneci6wHpjRJs+twE+95enAS13ddjC+etkmAuR4ywFgCXBSfx9Tf7ZJi/U/ABYBT/X38QyENgFKgYL+Po5Uvgb8mH5BQYEWFRX1dzWMMWZQWbZsWbkmuDlrwM+9U1RUhE3DYIwZTLbtr2VTWQ1nTB/db3UQkYQzGQz4oG+MMYPNx297hYirlN78mf6uShw7kWuMMR24/sk1FP/vi13Ku/blUnZc8zoRt+fD5u7+Pey45nV23di1MrvLevrGGNOBP79Z2uW8lYu3k9vDcj7cU8Xo3CDD9sSe+uMeyuzhnjpmPX1jjOnA4YXDOs/kyXUTPaataw7+5j3uueVN6OPnWSQt6IvIWSKyTkQ2isg1CdbPE5FK73rY5SJyXbLKNsaYvnLKlDzGDutaIM6LZvW4nIn4+FKDH+jbKyqTMrwjIj7gTuCTwA7gXRF5QlU/aJP1dVU9NxllGmNMKlzywV+5uj4AnJOS8mpD4T7df7J6+icAG1V1s6qGgAeB85K0b2OM6Tcj6gLURU9NuE7VZffux1CNdnl/taEIK3dUALBi++O8ue72NvvseV27IllBfzywvcX7HV5aWyd7t8Y/KyIz29uZiMwXkRIRKSkrK0tSFY0xpie8KByuj1uzufT3fLD2h6zbeGuX93bFovf53B1vUlUfpnzD1dTv/G3rSB8dHD39RANebb+v3iP2+K5ZwO+Ax9vbmaouVNViVS0uLIy7ocwYY1JMoKL1vU7RigpK190BwLZtf+nynlZvq2AcQkPYbU58/69Ni6qRXtW0M8kK+juAiS3eTwB2tcygqlWqWuMtPwMERKQgSeUbY0yfUA3iMhw35LZKr1u1Gqc21iv3Sdd751eFAjxMLm64xZDQ9iXN5fXxidxkBf13gWkiMkVEMoALgSdaZhCRsSKxa5FE5ASv7P1JKt8YY/pErfsJADbeFz/UHFzthdD9PmpefTVufaK5zWZHYttoi55+ZX1z716Cvapup5IS9DX2e+Ry4DlgLfCQqq4RkUtF5FIv2+eB1SKyArgduFAH+mxvxhjjya6KH8XOXBdLCy0bw/ZvX9qU/v+c1/iJ/wES3Zhb7sQSW67aVdF8vqAXN/N2SdLuyPWGbJ5pk3ZXi+U7gDuSVZ4xxvSnSFkZkVGx5YzT9rDrUy6Hvf8euWVzOMlZy6m+Vbiq+Nqc8nwr4FLU4KCB5vRIXXVzhj6O+nZHrjHG9EBo21Y04AXoYGyo5sCUpwF4zD2V28JfxE0QwH0NDQCs2XyQGnKoJI9Jh1Y2re/r4Q8L+sYYA1C2Hso3dm+bNiM+9SM2AfC2O5NH3dNo2FoVt8kFxKZ1GP7cO3xb7uMyubdH1e0pC/rGmAFvy/slPPO7XxIONcSt279oLWULVybYqmtUlS0/WUz1b2+EOz7SYd6td6/AbWi+6sZ3oOPpGfb9aRVVZbUJ1+WFapqWpUX/fvdr8V8UyWRB3xgz4B3YtYO1b7yCG2l9DXvDhm3Ur9xGw+ZKSh9+pUf7bqjYx32Zr3Kb/wge49Md5vVtqmLlS1ua3ksowf5ovnzTBXaXH0q4L0cTn1Ldsm8LH/h2dF7xHrKgb4wZ0GqrKnnl/j8BsO29Za3WZT5wLAWZF/H3zCeJvltL9f5yACr37cGNRtlXupnaqsrmDSq2Qcm91FZWcNsF57Jt9QqqDlY0rV7BDHjwPxLefdsoY8PzHdZ3j9O8v6tlBf+z6LmE+Ry3+VdCVfWfm5a3OPvY7RzssIzeSNrVOyJyFvBbwAf8SVVvbrNevPXnALXA11T1vWSVb4xJjqVLl7JkyRKuuOKKlJW596abyDsKgp/4GgwfR8h1EYTqPbt44CdXNeV74ne3wO9uoSCnljO/Np8JQKY0cDV/ZH3WKyy8LPZ4QgXUH8CJxHrducNzmffdyxh10wUER4XZl30DcAy3P/wQecOPa1WXLc+PJC/zl4z87I/An4ECYSdMhhsAIFh3AA2F2LV7CRXfiJ9zZ6dzANxx+Ijyz4xr+V3034CvgCra0HyVjk+br9OvGP8qI3aezm45SKVTSyW10EezMaRyls2zgWne60TgD97fPqGqSB/PS20Gnkgkgt/f8ce6vq6ezGBm0+cjHI3iRkNkZsSmxdVoFEQQp/mHsEZCiD8jvrxQA/6M+Idd7NtXxvCsLDKzMqhYMJl3OYlP3PAPREjK57Lt51s1dh+nI4LrujiOQyQcBp8PB2Lvd66g4p8/JX/ul5App8OwAkQcxHGIRiKEIxFEIjzx5jNE6n2Eoy4aOkRGVi71NQeoWP4IWUWnkllwONWuUhcJEfT5Kauv4+i8XMojSr0rjJQw2X4fTqgWskbEgt3rt6Gzv8Ib67bx9BsvEBkxhYLwAbJefZ5p4VxG7NjO3E+vhpU380u+xQPHfRq/L8Q5q95i+Wnnc2xVGR+tHcsxNYWsr1zBjZ+fyHWhSbScGOFIWcfVR68DYJtvNJOi+/j+tB/hW7Ob2Vkr2frwtbx29Gc45GTBrr2UzJrFiNEOQXc5hwqy8NULv9z9S945+hi+4fs+C6+/gWGs4M3p9fxm1JWU5kzk7efLKQ37OOXlVQRGf5NTVk/n8NxSzq9byqFcP+vlMC71X80X/TA/9H1KdQyfcd5h58+nsiujkC8e/2sunfA3fj1zPnBsU90rs3cwAqh0YuP/w7RvHqACIMm4P0pETgYWqOqnvfc/AVDVm1rk+SPwiqr+zXu/Dpinqrs72ndxcbF298HoqspVd/+Mb+14mGOc0m5ta4wxiRTVL2pa/k3gDhZHj+ef7qkcITt5KfM/u7SP992pnH3GPa3SAtrAf62/nY9uvZB8HcY+p5IcDZKvOUy4+WM9rq+ILFPV4rbpqZxls6szcfZ6ls1oOERuqNoCvjEmKZ6JntC0nEmIE50PGS+xWWQaCHR5Pw9FT49LC0smgrLJt4dMAkx0C8jXnN5Xuh3JGtPvyiybXckTS1RdCCyEWE+/u5XxZ2Tyv5f/gmj9DaiA+P2oasKf/aqKquI4TqvlpFHt2uPPuppvgOnuMFo0GsXn83WaL1xbi2RkdDpU09hurusiIoRCETIzA22yKA2ReqKu4hMfwYxM3HCY2v37yBk7noMVZQzPHYnP58N1m8dZHRQcH66rOE7sGA/sK6fk7Xc5/mPHM3rkuFZtoKpouB4nIwvXdYmGwkRCtWx77VHeOpDLlz5/PsGgH5zm43ejUcRxmtow0bBNXPuqQuV2IsMn4vfq5aoigLrR2P7dKI4v1nb1UZegL8Fnur6KquoqfDn5VIQjRJ0AIyTMbbf+Ggfl2ut+yqHqA2RkZJNZuZkPfKPIys5jfDCD+lCE2qpybj8YZXRmgG9OKCTbF7v3VIj9jy3hQ9SXbyMyfByZT1zG3llXMmLysWRnZbJr/YfsXbOEF597hUkVEbZPnswpOSsYtb+MsglTuOCE/0YRjv2whM98uJbMEROZECzishMn8I3tj/NK/kcI+QO8/e5X4g7r1vCXuET+yXWRr+Kvq2HinlJ2TJpOTk4uk7atpjKcwwuf+Bx78kYwc+dmqoPZTDqwl5A/gOvP5dLdo3A1yiGpI0t85I/Yy51Fx5OXUc9n33+L2R87k6yTP2Tfsz8k5BTwLf9NrD2Yw+7c4VQNz2DWKz/nefcBvu1/mpvDFzJtzkXMPW4cY6cdwc9cl29cfw9ZDdNxc97h7FOOZ78UUlR2CmOjh8UdS01VAznDkzvUMySHd4wxPbdt2zbWrVvHJz/5yZSVuf2zxdRsOMTkRYvInjObslAYQSjI8HPbBfEP28uta+DW71zHhiXN616pv4OSrY9w2gdbeH36JCaXV7K1IC+WX4Qjd+7luFnbqMgdzsFhR/NkZDYAC/77Gt66/lHecXYwpnoNn/vi16j8Z6zj4JdSIlrEXycqqwK7uWhrlBMv/yT54wp58+kfU5/1SFzdgi/eyhrnEeb7n+YR+Tyf/2nr4ZyNP72bYMN0DkzbwS1TSnmHuSx+XpHGgRffKojGxvsP+/ncHndC2xveSVZPv2mWTWAnsVk2v9wmzxPA5SLyILETuJWdBXxjTOpNmjSJSZMmpbTMiU+27tgVZjT/Wptz9ud479lWk/bykdI9rD3pGLYUPkvFI1vIZhSj5sKVl/8eJzuLosqDjCwYzW+/+SUmzjyOL173cwBWbNnDbxZv4eq548i97nL+7T/OgUCQk2/8MtOun0IgL0Luyffz9/2vc8Tbf2Hq5b9GfruCr2wXINYTF4kF4aOy57Ly0CNodnO9ArWjWeLfwBStpjOC8v94iLN5moPOlYx0Y78cfURpvCYoqaMOnqQEfVWNiEjjLJs+4N7GWTa99XcRm4ztHGAjsUs2v56Mso0xQ9sZX5vP9FNPp2LvHo6e23pMfErxKexb4setVcZ+tvliwNEjRwLwnbsfICPY/LDyWVPG8ucpYwE45u9PNaWLCIU0Xxv/zXM/BufGTqLulj1EdWxzXn9sSC2TXMb9MINdv4/doeU/lE9WxTSqnDr2RGMzsa3j6ARH1Dxcl0cVeVQ1pfhlHc6ZhxN9HhQ3wba9l8pZNhX4brLKM8akj3FTj2Lc1KMSrhv93RMSpgNkD8/reiHn/gbGz4lLDsiW1kHf631rm7uDVaKEvUDdoLHLe+sli7aaBtQ1PqgLIfIKg5QBIvFTTiRD0oK+McYMasXtDT60PpHueMM7RFvfmKVuPaqxvJXu4dwTOZuGjERPRHG8vTafT61gOPlAOcMpyMgCoojTNxd22DQMxhjTgXrJbvXe8TX29JuDfvZrDmiI6L4NAJS507kxchG1Tvyll5pgqYHYL4NazcSZehLB0VWM+nLiXza9ZUHfGGM6EJHWd2KLd8mtelM8EAGnFtSBQE0l/1H/MXKJfVEk7qvHp0qLNPE5FPzgMwRn9s3JdAv6xhjTgf3+Ma0TnNiouHj3QTTdmOCARCJkkUFjYHcS3MMSCuwDQP3Nk+s0DvWk4vmxFvSNMaYDlb5Rrd43DulnTPZ64o1Bn+a/jadoE/X0wz5vvnx/GPVCsHo5o+38NkimXgd9ERkpIi+IyAbvb347+UpFZJWILBcRu9vKGDMoNQ3FeL34F1cH+K/ZGa0ivDbljact7r6ee/IrzJn9QFPw382oBFskVzJ6+tcAL6nqNOAl7317zlDV4xPdJWaMMQPR1qzWobtx6gvJzCTjiCPYVO3n32Z8J7ZSW/1JGPQ3Z80CoCL/GLKyxpOff1IK+vfNknHJ5nnAPG/5PuAV4MdJ2K8xxvS7V0f5mV4eYiKxE7iOP9ZXDh55JEc8/RQPefle3v3rpm2agn6CbvWO4UewdHcl2fmHx60bLGP6YxqnU/D+jm4nnwLPi8gyEZnf0Q57O8umMcYky5SCYWR4ffFdk4a1f/18U7JLtRe+6934vF+/4Fg2zBvHJ06a2JQW9rKtJhKXP9m6FPRF5EURWZ3gdV43ypqrqnOIPUzluyJyWnsZVXWhqharanFhYYGbQqEAABO+SURBVGE3ijDGmOT66WdnkhuI9fLdDiKmZkJkTGzyhOe8x169dij+5qy87AA/Pms6/hYznzZex7O2j6ZeaKlLwzuq+on21onIXhEZp6q7RWQcsK+dfezy/u4TkceAE4DXelBnY4xJqexwrOc+NjfRHbbNGo5V3N1CNvUUSCV1MqJL+1+TDbOroT6n7ydJSMbwzhPAxd7yxcA/22YQkWEiktu4DHwKWJ2Eso0xps9lHR8bcTjs9Imd5Ixd3XOBs5LXMr/P+cO2dWn/7+T5OIsqhh+W26t6dkUyvlZuBh4SkUuAbcAXAETkMGIPSD8HGAM85j0Mwg8sUtV/JaFsY4zpcyO/eBSRj08iMDq73Tw+ZzgTD/s6PoF/d4/h1w1f5rJLPtul/asINbR3B29y9Troq+p+4MwE6buITaWMqm4GZvW2LGOM6Q/iSIcBH2DevPcBqPrUNngOrvrZ77s8aVpuZiwUf+a4cb2raBfYLJvGGJNEw8+YxPAzujdvjs/7cijMTe6jEROxaRiMMSaNWNA3xph+9u3TYjdqzZrQtat9esOGd4wxpp+dMrWA0ps/k5KyJPYUw4FLRMqArf1djwGkACjv70oMMNYmrVl7xEvHNpmsqnF3tw74oG9aE5ESm7CuNWuT1qw94lmbNLMxfWOMSSMW9I0xJo1Y0B98FvZ3BQYga5PWrD3iWZt4bEzfGGPSiPX0jTEmjVjQN8aYNGJBfwAQkYkislhE1orIGhH5XoI8+SLymIisFJGlInJMV7cdjHrZJkHv/Qpv2+tTfwTJ15s2abHeJyLvi8hTqat53+ltm4hIqYisEpHlIlKS2tr3E1VN2Qu4l9hDVlanstyB/gLGAXO85VxgPTCjTZ5bgZ96y9OJPYy+S9sOxlcv20SAHG85ACwBTurvY+rPNmmx/gfAIuCp/j6egdAmQClQ0N/HkcpXSk/keo9IrAHuV9VjOssPUFBQoEVFRX1aL2OMGWqWLVtWrgnuyE3p3Duq+pqIFHVnm6KiIkpK0uNXlzGd0Zf/l/sb8nl0zKdYWnmIVXNnUh6K8Mcnf81V+VE+5KMElt1NJCePHxZdREY0zEuH1zDi8a93uYyVOdM4rmZDt+sWFh8BjXZ7u+56ceRJfOLAOwnX/WrSV/nBtvu7tJ9vzVjA3R8sSGLNOnbRzJ/zQsHcuPTCqoP8+/uvJtxmwYIFPS5PRBJOXzMgJ1wTkfnAfIBJk7o3L7UxQ5m8disXAz8+PRY8Fm4v4+UDVby07hcA/IXvs4AX4SBcfNyVAIx4/AvdKqMnAR9IScAH2g34QJcDPpDSgA9w4d5nEwb9suH5Ka3HgDyRq6oLVbVYVYsLC+N+nRhjPHaXzeDxXu6M/q4CMECDvjHGDDUqqXgCbucs6BsziFlP33RXSoO+iPwNeBs4SkR2iMglqSzfmKFIGBg9SNMxHSD/Tqm+eudLqSzPmKEuhVdcm14aKEHfhneMGcR0wIQSM1hY0DfGmBQYKD/KLOgbM4gpWE9/kLCrd4wxxqScBX1jBrGBMmRgOjdQzr5Y0DdmsBsYscQMEhb0jRnMrKs/aFhP3xhj0shA+X62oG+MMWnEgr4xg5g2/ccMeHbJpjHGpA8b0zfG9JpaN3/QsKBvjEmOgRFLzCBhQd+YQcz6+YPHQPm3sqBvjDEpYHPvGGN6zebTN91lQd+YQcxivukuC/rGGJMCdvWOMabXrKc/eFjQN8YkxcAIJWawsKBvzCBmPX3TXRb0jTEmBVy7ZNMYY0yqWdA3ZhBTtdl3TPdY0DdmkBsYgwZmsLCgb8wgZr18010W9I0xJo1Y0DfGmDRiQd8YY9KIBX1jBjEb0zfdZUHfmEHMplY23WVB35hBTuyiTdMNFvSNGcTs1izTXRb0jTEmjaQ86IvIWSKyTkQ2isg1qS7fmKHE+vmmu1Ia9EXEB9wJnA3MAL4kIjNSWQdjjEln/hSXdwKwUVU3A4jIg8B5wAfJLuilG89gfLQs2bs1pl9N9/4ufvWihOsv5JGm5fbymP6h5SEY171t6uvDBIOBpNYj1UF/PLC9xfsdwIltM4nIfGA+wKRJk3pUULlvxIB5PJkxyTI9up0agpT6Y9HDrxEUYbRbwSitYoevkKLoHoKEm/JkRkMcrnv6s9opU00WudR1Ke86ZyJHuds7z5gk/8g/M2H6/L1/Agri0mdFCvvkyizRFF7oKyJfAD6tqt/03l8EnKCqV7S3TXFxsZaUlKSqisYYMySIyDJVLW6bnuqe/g5gYov3E4BdHW2wbNmychHZ2sPyCoDyHm47VFmbxLM2iWdtEm+wtcnkRImp7un7gfXAmcBO4F3gy6q6po/KK0n0TZfOrE3iWZvEszaJN1TaJKU9fVWNiMjlwHOAD7i3rwK+McaYeKke3kFVnwGeSXW5xhhjhv4duQv7uwIDkLVJPGuTeNYm8YZEm6R0TN8YY0z/Guo9fWOMMS1Y0DfGmDQyJIN+uk3qJiKlIrJKRJaLSImXNlJEXhCRDd7f/Bb5f+K1zToR+XSL9I94+9koIreLyKC5pVlE7hWRfSKyukVa0tpARDJF5O9e+hIRKUrl8fVEO22yQER2ep+V5SJyTot1Q7pNRGSiiCwWkbUiskZEvuelp9fnRFWT9gJKgVXAcqAkwXoBbgc2AiuBOcks3yvDB2wCDgcygBXAjGSXM5BeXrsXtEm7BbjGW74G+IW3PMNrk0xgitdWPm/dUuBk79/pWeDs/j62brTBacAcYHVftAFwGXCXt3wh8Pf+PuYetskC4IcJ8g75NiE2880cbzmX2D1DM9Ltc5LUE7kiUgoUq2rCu9a8XsUVwDnE5tz5rarGzb3TUkFBgRYVFSWtjsYYkw6WLVtWrqqFbdNTfZ3+ecD9GvumeUdERojIOFXd3d4GRUVF2Nw7xpihbOWL7/KPN57msq98i9FTxydln+1NX5PsMX0FnheRZd5MmW0lmmUz7ghFZL6IlIhISVmZTY9sjBnaVq1YBcD2taV9Xlayg/5cVZ1D7CEp3xWR09qsT3RiMG58SVUXqmqxqhYXFsb9OjHGGNNDSQ36qrrL+7sPeIzYQ1Na6vYsm8YYY5InaUFfRIaJSG7jMvApYHWbbE8AX5WYk4DKjsbzjTEmnSTzwpr2JPNE7hjgMe9yVT+wSFX/JSKXAqjqXcQmWjuH2CWbtcDXk1i+McaYTiQt6GvsubezEqTf1WJZge8mq0xjjBlKUnE/5JC8I9cYYwajVAzvWNA3xpj+lsIJTyzoG2NMf0vhDPcW9I0xA1rNwSpKV27o1zpEI2E2rnwrLr2s7iAfHNjUDzXqOQv6xpgB7d7v/4BHf/b9fq1DyZ9/yNR/nM2GVUtapX/8nff5+Irq3hdgwzvGGBMTrtvT31UgZ/9KAGrKd7RKL2Nkcgqw4R1jjBk4NJVd8T5mQd8YY7qozy6ptOEdY4wZgFJwHX1fs6BvjDGdabxTVt3+rUcSWNA3xphONI7pD/5+vgV9Y4xJKxb0jTGDQirmpWmf19MfAl19C/rGmMFhAARcwcb0jTEmJfoy5u/dtoGlj/4mLv3Dg1u4ZfEi3j70iVgd2unqv7N3NX/f9AqbSz5kzeL3+rCmvZfMh6gYY0zfUaWvLmiP/vlcTtA9VJ55MXkj8pvSz16+mzpmcKmznjL3yHa3P/+DCDCCS199EICZZ8zpk3omg/X0jTGDgtuHA+rDtRKAqNt6+KaObG9JiRJAhsCgflKCvohMFJHFIrJWRNaIyPcS5JknIpUistx7XZeMso0xacLty4Drnah1o+2sV1xkIJxW6LVkDe9EgKtV9T3v4ejLROQFVf2gTb7XVfXcJJVpjEkjfXn1TuOeOy5DhsTlO0np6avqblV9z1uuBtYC45Oxb2OMAfr0TG7TzVft3XErQ2fStaSP6YtIETAbWJJg9ckiskJEnhWRmcku2xgzdLl9OAVCU0Bv94sldhJZ+3qAJwU/JJJ69Y6I5ACPAlepalWb1e8Bk1W1RkTOAR4HprWzn/nAfIBJkyYls4rGmMEqBUMr7QX1WKoN77QiIgFiAf8BVf1H2/WqWqWqNd7yM0BARAoS7UtVF6pqsaoWFxYWJquKxphBrC/P4zYP77RXSOwkbl/19CVuoe8k6+odAe4B1qrqr9rJM9bLh4ic4JW9PxnlG2PSQB9G/ebhnXau3nH67h6BWPltF/pOsoZ35gIXAatEZLmX9l/AJABVvQv4PPAdEYkAdcCF2r+TaRhjBpF2T7ImY9+Nf9stQmJfDEMgZCXr6p03VFVU9ThVPd57PaOqd3kBH1W9Q1VnquosVT1JVeMfLW+MMe3423W/4K1Fv+bpL3ysVfqaxe9x781/wI12/qVwcN9ONt44h12l69qsaT118luPvMxDv/m/prWK8hc5A+1CGZ358I2V3P3zO4mGm39VtBze2bVuG3fe+BvqKg71uqxE7I5cY8ygULH7PfJvWMjhq8pbpT/6ypNsq99LtCHS6T42vngPU6ObKH36tlbpbS/ZfH71a3xQsalpgjVxICoutHfvVjc89uKT7AyVUXegOn6lwktPPk9ZtIIP31zZ+8ISsKBvjBl0Wo4MS9PdtF3ohYsX8tqM43R2Irf55q3u1TPxvrxrgZzE4dc79dln005Y0DfGDDoRt0WvvukcbHeCfjsBtU163MQLSQn6nnaib7e+xHrAgr4xZtBpGRAbg6Tbhat7mnrX7fT026Y3Bn2VxkCcjN53Zz19L5f19I0xJsZt0dNvPAnapSDp9fTbPgylKai3c/mOiBf8k9jTdyTxJaCS1C+YeBb0jTGDTjQSblpuGg7pwvCOtDOm36htUG/q6Td+tSQhEDfd4JUg6CvN9wNYT98YYzyu2/Jyx673jLWzE7l0PKafjEDc2T4cx4K+McYANA3KuNH4aye7cvNWY0+/7cNQOuvJN34ZJHN4J1FQj32BWdA3xhgAtDE2txrT7/oljs0nTxMH/bY9fafppGvT2dXuVrlbFG0a0++rsobsM3KfePjPhBrq+rsaxpheio5tnpdxbSALgPcfW0RGZuxRhvWZNahE+dczD+DL9HW4r5q9W9nJqeyvzmTHX3/flC46kwyOpOLFJ8h6ezEajMWOvOrYYxRdfwgnqKzfUcqeFtuNHj611f4bt3ukRZ5GTvYwhudl05AT+5XyyptPEQg6BH1B6hrqAajcX05dfS0AFQcOoNriSyBJZKBPf1NcXKwlJSXd3m79ghkcyc4+qJExxnTf9Vlf5w8nfC0u/Yy6lzlqaduZ6GG4m8Vl1/6QYKDjL7L2iMgyVS1umz5ke/r/yj6NFwn1dzWMMb2k3n+aOrxxE15K0yNOkie2N1eE+kCAzHA9PnXiyqjMzKbayWFsQzkhn4/scKTdutRl5fDlXf9CcYiKg19jQ1RjD9VwZM5YKuuijMj2gSrltVEmFg7D7yS3lw9DOOhf+aO7+rsKxhgz4Az44R0RKQO29nDzAqC801zpxdoknrVJPGuTeIOtTSaratxTqAZ80O8NESlJNKaVzqxN4lmbxLM2iTdU2sQu2TTGmDRiQd8YY9LIUA/6C/u7AgOQtUk8a5N41ibxhkSbDOkxfWOMMa0N9Z6+McaYFizoG2NMGhmSQV9EzhKRdSKyUUSu6e/69DURKRWRVSKyXERKvLSRIvKCiGzw/ua3yP8Tr23WicinW6R/xNvPRhG5XZI96UcfEpF7RWSfiKxukZa0NhCRTBH5u5e+RESKUnl8PdFOmywQkZ3eZ2W5iJzTYl06tMlEEVksImtFZI2IfM9LT5/PiqoOqRfgAzYBhwMZwApgRn/Xq4+PuRQoaJN2C3CNt3wN8AtveYbXJpnAFK+tfN66pcDJxG5yfxY4u7+PrRttcBowB1jdF20AXAbc5S1fCPy9v4+5h22yAPhhgrzp0ibjgDneci6w3jv2tPmsDMWe/gnARlXdrKoh4EHgvH6uU384D7jPW74POL9F+oOq2qCqW4CNwAkiMg4Yrqpva+zTen+LbQY8VX0NONAmOZlt0HJfjwBnDvRfQu20SXvSpU12q+p73nI1sBYYTxp9VoZi0B8PbG/xfoeXNpQp8LyILBOR+V7aGFXdDbEPOjDaS2+vfcZ7y23TB7NktkHTNqoaASqBUX1W8751uYis9IZ/Gocx0q5NvGGX2cAS0uizMhSDfqJv1KF+XepcVZ0DnA18V0RO6yBve+2TTu3WkzYYKu3zB+AI4HhgN3Cbl55WbSIiOcCjwFWqGj+vcYusCdIGdbsMxaC/A5jY4v0EYFc/1SUlVHWX93cf8BixIa693k9QvL/7vOzttc8Ob7lt+mCWzDZo2kZE/EAeXR86GTBUda+qRjX2bMG7iX1WII3aREQCxAL+A6r6Dy85bT4rQzHovwtME5EpIpJB7ETKE/1cpz4jIsNEJLdxGfgUsJrYMV/sZbsY+Ke3/ARwoXeFwRRgGrDU+0lbLSIneeOPX22xzWCVzDZoua/PAy97Y7mDSmNg8/wbsc8KpEmbeMdwD7BWVX/VYlX6fFb6+0xyX7yAc4idld8E/Hd/16ePj/VwYlcXrADWNB4vsTHEl4AN3t+RLbb5b69t1tHiCh2gmFgQ2ATcgXfH9mB4AX8jNlwRJtbTuiSZbQAEgYeJnchbChze38fcwzb5P2AVsJJYcBqXZm1yKrGhlpXAcu91Tjp9VmwaBmOMSSNDcXjHGGNMOyzoG2NMGrGgb4wxacSCvjHGpBEL+sYYk0Ys6BtjTBqxoG+MMWnk/wNY1D750VLQHwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "cycles = [] # names of cycles with data\n",
    "for filename, Di in D_2l.items():\n",
    "    cycles += [str(Di['cycle']).zfill(2)]\n",
    "cycles.sort()\n",
    "    \n",
    "# x1 = 2.93e7\n",
    "# x2 = 2.935e7\n",
    "\n",
    "beams = ['gt1l','gt1r','gt2l','gt2r','gt3l','gt3r']\n",
    "\n",
    "### extract and plot data from all available cycles\n",
    "fig, axs = plt.subplots(4,1)\n",
    "x_atc = {}\n",
    "h_li = {}\n",
    "h_li_diff = {}\n",
    "times = {}\n",
    "\n",
    "for cycle in cycles:\n",
    "    # find Di that matches cycle:\n",
    "    Di = {}\n",
    "    x_atc[cycle] = {}\n",
    "    h_li[cycle] = {}\n",
    "    h_li_diff[cycle] = {}\n",
    "    times[cycle] = {}\n",
    "\n",
    "    filenames = glob.glob(os.path.join(datapath, f'*ATL06_*_{rgt}{cycle}*_003*.h5'))\n",
    "    for filename in filenames:\n",
    "        try:\n",
    "            for beam in beams:\n",
    "                Di[filename]=atl06_to_dict(filename,'/'+ beam, index=None, epsg=3031)\n",
    "\n",
    "                times[cycle][beam] = Di[filename]['data_start_utc']\n",
    "                \n",
    "                # extract h_li and x_atc for that section\n",
    "                x_atc_tmp = Di[filename]['x_atc']\n",
    "                h_li_tmp = Di[filename]['h_li']#[ixs]\n",
    "                \n",
    "                # segment ids:\n",
    "                seg_ids = Di[filename]['segment_id']\n",
    "#                 print(len(seg_ids), len(x_atc_tmp))\n",
    "                \n",
    "                # make a monotonically increasing x vector\n",
    "                # assumes dx = 20 exactly, so be carefull referencing back\n",
    "                ind = seg_ids - np.nanmin(seg_ids) # indices starting at zero, using the segment_id field, so any skipped segment will be kept in correct location\n",
    "                x_full = np.arange(np.max(ind)+1) * 20 + x_atc_tmp[0]\n",
    "                h_full = np.zeros(np.max(ind)+1) + np.NaN\n",
    "                h_full[ind] = h_li_tmp\n",
    "                \n",
    "                \n",
    "                x_atc[cycle][beam] = x_full\n",
    "                h_li[cycle][beam] = h_full\n",
    "                  \n",
    "                                            \n",
    "#                 ### here is where you would put a filter\n",
    "#                 # you would want to de-mean and detrend that section first:\n",
    "#                 h = h_full\n",
    "#                 x = x_full\n",
    "#                 h = h - np.nanmean(h) # de-mean\n",
    "#                 h = scipy.signal.detrend(h, type = 'linear') # de-trend; need to deal with nans first\n",
    "#                 # use scipy.signal.filter to filter\n",
    "\n",
    "#                 # differentiate that section of data\n",
    "                h_diff = (h_full[1:] - h_full[0:-1]) / (x_full[1:] - x_full[0:-1])\n",
    "                h_li_diff[cycle][beam] = h_diff\n",
    "\n",
    "                # plot\n",
    "                axs[0].plot(x_full, h_full)\n",
    "                axs[1].plot(x_full[1:], h_diff)\n",
    "#                 axs[2].plot(x_atc_tmp[1:] - x_atc_tmp[:-1])\n",
    "                axs[2].plot(np.isnan(h_full))\n",
    "                axs[3].plot(seg_ids[1:]- seg_ids[:-1])\n",
    "\n",
    "\n",
    "\n",
    "        except:\n",
    "            print(f'filename={filename}, exception={e}')\n",
    "\n",
    "\n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "oops! Not all cycles in beam gt1r are the same length. Trimming\n",
      "oops! Not all cycles in beam gt3l are the same length. Trimming\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEFCAYAAAAPCDf9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeZwdZZ3o/8+3Tp2lz+k93dmXDiEBww4tLnEUZFDABZdRYRzEkTt4XcaZOzP3CsyM46iMXH/KjI4jI6jjMoLiVQavIAgIFxEQEhbDkpBA9nTS6aT3s1d9f39UndOnl3SW053udH/fr1e/zqmnnqp6TiX9PU8/9SyiqhhjjJkdnKkugDHGmGPHgr4xxswiFvSNMWYWsaBvjDGziAV9Y4yZRdypLsChtLS0aFtb21QXwxhjjivr1q3rUtXWkenTPui3tbWxdu3aqS6GMcYcMx//4ToKnvK1y86iJhY5qnOIyLax0qd90DfGmNlEVbl7/R4A0vniUQf9gzlkm76IfEdEOkXkuYq0ZhG5T0Q2ha9NFfuuFZHNIrJRRN5akX6OiKwP931NRGRCP4kxxswAX7jrxfL7ObXxCT//4TzI/S5w0Yi0a4AHVHUl8EC4jYisBi4DTgmP+YaIlL6mbgKuBlaGPyPPaYwxs5LvK6/sG0BV+fYjWwC46JT5k3KtQzbvqOrDItI2IvlS4Lzw/feAh4BPh+k/UtUcsEVENgPnishWoF5VHwMQke8D7wJ+WfUnMMaY49zdz3XwyVufHpb2sfNWTMq1jrZNf56qdgCoaoeIzA3TFwGPV+TbGaYVwvcj08ckIlcT/FXA0qVLj7KIxhhzfOhJF8rvP3/pKVx+7lLcyOT0qJ/os47VTq/jpI9JVW9W1XZVbW9tHdXjyBhjZgxVpWsgB8D8+gRXvK5t0gI+HH1Nf6+ILAhr+QuAzjB9J7CkIt9iYHeYvniMdGOMmbWe2dHDNT/9PRv29ANw4wfOmPRrHu3Xyc+BK8P3VwJ3VqRfJiJxEVlO8MD2ibApqF9EXhv22vlQxTHGGDMj7OxOc+3Pfk+24B00j+8rHb0Zfvzkdt5706P0pAu89ZR5fPYdq3n9ipZJL+Mha/oichvBQ9sWEdkJ/ANwA3C7iFwFbAfeB6Cqz4vI7cALQBH4hKqWPv3HCHoC1RA8wLWHuMaY48I3HtpMXSLKFa9dNm6+f/9/L3PbEzs4Z1kzf3TO4jHzvOoz95Ar+gCsOXEON/3JOdQnohNe5oM5nN47lx9k1wUHyX89cP0Y6WuBU4+odMYYM8VUlS/dsxHgkEG/Nh4E7909mTH3r9/ZWw74J7SkuOVD7SRjx3aMrE24Zowx4+iu6FnTmymMkxPibhBS73th76h9dzy9k3d947fl7V986g3HPOCDTcNgjDHjemRzV/n9/oEcDTUHb4rp7M8CsH5XLwC/2bSPf3twMzsOZNjVk+F1J8zha5efRWvdxI+0PVwW9I0xZhy3/W57+b1/iCXFb3tiR8X77Vz7s/WkYhEueNU83n3WIv78ghOJuxM7l86RsqBvjDEhVeXBjZ185LvBzL5bb3gbmzr7y/sPDObHPf6tp8zj3ueDpp1rf7aepc1JfvkXf0AqPn1CbVUlCadX6Ac8oKiq7SLSDPwYaAO2Au9X1e4w/7XAVWH+T6nqvdVc3xhjJorvK//t+2v59YbOctpND71M18BQoP+3Bzezce88ntvZi6fKeSe1ctEp88uDqdbv7KW1Ls7bTlvA2cuaeMOJLdMq4AOI6iH+Xhnv4CDot6tqV0Xal4ADqnqDiFwDNKnqp8PJ2G4DzgUWAvcDqyq6dI6pvb1dbT59Y8xk2tI1yEd/sJaX9g7wFxes5MLV83j7vz4yLE9zKlau6TcloyjB9Akx1+H1K+bw0Teu4PJbHueiU+bz71ecMwWfYjgRWaeq7SPTJ+Mr6IgmYwMem4QyGGPMYfF85fwvPwTAa5Y38xcXrMRxhF/8+Rt4ZkcPqoobcdjaNcg3H36F//jwqznvpFZ8hR89uZ2/veM5Htq4j4c27gPg3WcfdFqxaaHaoK/Ar0REgW+q6s0c+WRso9iEa8aYiaKq7OvP0ZiMsX5XD++96TFOnl/HH75qHp9884lEnKGpwb73kXNxwu1TFzVw6qKG8r5sweMdZywsp0UEPviaZbxpVStrt3bzL/e/xNb9aV63Ys6x/YBHqNqgv0ZVd4eB/T4R2TBO3sOedC388rgZguadKstojJllPF954MW93Pnsbh7Z1EVvpkA0IhS8IJxs2NPPhj39eKq0hAuV/MeHX00ievCeNYloZNiXQMnipiSLm5K866zpXcMvqSroq+ru8LVTRO4gaK450snYjDHmoIqez87uDPsH8xQ8n9a6OEuaksTCgVCqymDeI+E6RBxhIFfkoz9Yx6Mv76elNsZbT5nH6gX17OnLUfB83riqlTetauVj/7mOmx56GUfgzCWNvGnV7JjR96iDvoikAEdV+8P3bwE+x9BkbDcwejK2W0XkRoIHuSuBJ6oouzFmBtvaNci//nozd63fTbbgD9tXl3C57NVLaKmN88VfDjUwxFwHR6DoKde/+1Q+0L7koNMUX3fJqyh4ypLmGv7HhavKzTozXTU1/XnAHeFSty5wq6reIyJPcuSTsRljZgnfV57d2cP6Xb30ZQq01sVZ1Jjk9CUN1Cei7DiQ5l9/vYmfPrWLaER4z9mLOWtJI611caIRhz29We57YS/femQLlZ0P33b6AubXJ8gUPN591iJe3dY8bjmWNCf51pWjOrfMeFV12TwWrMumMdNXrujRky6wty/L7p4sO7vT7DiQZldPhlzRxxEh4giOCHNSMUTgN5u62DXGhGQRRzixtZaX9w3gOMKfvGYZ//28E5hblxjz2tmCRybv0ZiMElY+TYVj2WXTGHMcyOQ9+rMFEFCFdN5jMFekP1tkMFdkMF8kk/cYyBXpzQSBvbM/x77+HF0DOboHC+Q9f9R56xIui5uSJKIOvq/4CgXP59mdPRQ9nzOXNPLXb1nFmhNbaKiJsq8/x/YDaR5/ZT/P7OjhvJNa+cgbljOvfuxgX5KIRsZ98GrGZkHfmBkuX/Tp6M2wuXOADXv6+f3OHp7b1TdmbftgHIE5tXHm1ceZWxdn9YJ6mmtj1Cei1NdEmV+fYEFDgiVNSRqSRzY3/JLmJEuak6w5cfIXEDEW9I2ZllQVz1e8sPnV85VswScd1r4H8x7pfJF0ziNd8BjIFtk/ENTAuwbydA3k6EkXOJAO3le24rbNSXLOsiYuP3cJjckYShDUk7EIqZhLbdwlFf4kYxFS8SAtMksedM50FvTNIakqqiBCue00X/TJFDx8X1HAV8UP86mOsx2ezw/T9CCvIkJdwqU+EaU27uJGBNcRVKHoBwGx6Pv4PhR9P9zWilefoq8UPR2xzydb8Ojqz3PKonpWzq0j7/lk8kH7cMH38X1lIFdERIhGhHTe46lt3Ty5tZtt+wdZNifFCa0p+rNF+jIFejJ5ejMFejMFip5SE4sEgTPmMqc2CKq5gke24FPwfGKuQ9yNkIg6xFyHXNGnL1OgL1ukP1OgL1ugL1Mcs+nkUOoTLi21cebUxmhrSXJWspF59QkWNdWworWWVfNqqTuGqzSZ6WfGBv2/+cmz7OxOB0EHQEHR8rYfBh4Ng5HnB6PHRBgdiBi+7avi+8ODl6/g+X55tFllzcp1godZpbZNDc9ZylwKhEPvS+lDJxEET5WRD95FBEfAERkKsqWdY5yndK7wTTlf6d74FecXEbwx5pIt3aPZZnlLihWttbyyb4CHX9pHfY1LQ02Uhpooc+sSnNhaixtxyIQ178Fcked39yECiTDIuxGHgVyRroE8uYJHrugTjzrUJ4LzLGmqob4mSn0iSioWKXcjjDhCwnVIxlxqYhFS8QjJWFATT8Yi1MajNKWiUz5tr5n+ZmzQ9/0gMCNhMHdAcMLaahAknYqAWarBqmp5f2W+ofdD+Z2KfCJBTbTyL2ARQVUp+KVgLcQiQ9cSGQrAIkMxuLJGLQx9KThhL4iKWD30ReQH+6Ucz4fOWzpP6RgYHbSDzxccVwrqihIRKQeeoS9QJeY6JKKRcs8MJ/zGdCruU+lcjgiOM2K7fE9LX1xBiR0n2C7VtnszBQZzHkXPp+AH5XEjQf7Sl6kbCV8dIeI4Q+nha6S8HQzeibkOjcko63f2sv1AmrjrkIxFqIm5RMNzp+JB8MwXFdcRzgi7DBpzvJv2XTZFZB+wbarLMY20AF2HzDW72D0Zzu7HaLPxnixT1VHDjKd90DfDicjasfrezmZ2T4az+zGa3ZMhtjC6McbMIhb0jTFmFrGgf/y5eaoLMA3ZPRnO7sdodk9C1qZvjDGziNX0jTFmFpn2/fRbWlq0ra1tqothjDHHzO7eDPsH8py6sKE81uZIrVu3rmusLpvTPui3tbVhUysbY2aT0z57L7FskZ/+zXm0taSO6hwiMub4pkM274jId0SkU0Seq0hrFpH7RGRT+NpUse9aEdksIhtF5K0V6eeIyPpw39fEJsA2xphR7nluD/3ZIm9a1XrUAX88h9Om/13gohFp1wAPqOpK4IFwGxFZDVwGnBIe8w0RKU0GchNwNcEyiSvHOKcxxsx6//0/1wFw3kmTs2bvIYO+qj4MHBiRfCnwvfD994B3VaT/SFVzqroF2AycGy6QXq+qj2nQXej7FccYY8ys5/lK2zV3lbf/dM3ySbnO0fbemaeqHQDh69wwfRGwoyLfzjBtUfh+ZPqYRORqEVkrImv37dt3lEU0xpjjQ086z4rr7i5vf+btqyftWhP9IHesdnodJ31Mqnoz4WCK9vZ2G0hgjJnR9vblyu8fv/YC5jeMv1RkNY426O8VkQWq2hE23XSG6TuBJRX5FgO7w/TFY6QbY8yslSt6nPR395S3v3b5WZMa8OHom3d+DlwZvr8SuLMi/TIRiYvIcoIHtk+ETUD9IvLasNfOhyqOMcaYGaHg+Xzxly/y4MbOQ+Z99zd+OyzgA7xp1eQ8vK10yJq+iNwGnAe0iMhO4B+AG4DbReQqYDvwPgBVfV5EbgdeAIrAJ1TVC0/1MYKeQDXAL8MfY4yZMVxH+N6jW/E85fyT5o6ZR1W59/m9PL29p5x2x8dfT0ttnIaayV/KctrPvdPe3q42OMsYM5W+dM8GohGH/3HhqkPmfW5XLwsaEsypHb3S2p7eLK/94gPD0rbe8LYJK2clEVk31hoC035ErjHGTLUtXYOk4ocXLk9d1DBm+s+e2slf3f5sefvHV7+W15wwZ0LKdyQs6BtjzCHc9Cfn4PuH1ypy5zO72LY/zacuWMl//HYL//h/XxiVp31Z05QEfLCgb4wxh8VxDm/mmP/7bAf3v7iXG+97acz9k9Wcc7gs6BtjTIXSc87K6cF60nnuXr+HNSfOYdmc8efDeezl0euv/+jq1/KqBfXH5EHtoVQV9EVkK9APeEBRVdtFpBn4MdAGbAXer6rdYf5rgavC/J9S1Xurub4xxkykdL7I27/2CK90DQLw22vezJobfj3uMZuuv5hoZKj3+//+o9P55K1Pl7enumY/UlW9d8Kg366qXRVpXwIOqOoNInIN0KSqnw4nY7sNOBdYCNwPrKro0jkm671jjDlWKue+ORLLW1L87/eezrnLmwHwfT3s5qDJcrDeO5OxctYRTcY2Cdc3xpgjsqVrkP/v3g0ArF5Qz6brLx62/18+cOYhj3//Nx+j7Zq7+PK9G6c84I+n2jZ9BX4lIgp8M5wzZ9hkbCJSORnb4xXHjjvpmjHGTBTPVwRIFzxO/YehVuU7Pv56zlraxPlffqicdten3oCIsOWLlzCY96gNu2q+66xFeL4SGRHQL/jKQ5y1tIn/sy6YU/LrD27mb9560qR/pqNVbdBfo6q7w8B+n4hsGCfvYU+6JiJXE8y9z9KlS6ssojFmNvrCL17gW49sGTfPu7/xKFu+eEl5e+MXLio/wBWRcsAvGRnwAR746/MA+NM1bbzta4/w9T8+q8qST66qmndUdXf42gncQdBcszechI3DnIxtrPPerKrtqtre2jr5c1EYY6a3fNFn/0CO53b1oqqMfBY5crvtmrvGDfgfO29F+f3ya4emNI67kbGyH5ZTFjaw9Ya38fbTFx71OY6Fo67pi0gKcFS1P3z/FuBzDE3GdgOjJ2O7VURuJHiQuxJ4ooqyG2NmsCe3HuB9//5Y1ee5/6/eyIlz68gXfXxVEtEgsN/88Ct4FQOu7vzEmqqvdTyopnlnHnBH+KeQC9yqqveIyJMc+WRsxphZZDBX5JKv/YZt+9MAfOqClfxVOK+N7ysnVCwocrT+6xNrOHNJY3k75g5v2Hj5ny7h+d29tNTGmVc/udMZTyc24ZoxZkJ09mX56gOb2NWT4bxVrfRkCtz6u+109uc4a2kjCxoS3L1+zxGfd+3f/SEtIyYve2lvP2/554fL2x87bwUffeMJ1CeiFHyfWMQZNrhqNjpYl00L+saYsqC9HAbyRYqe8p+Pb+OrD2wa1gwylrjrkCv6E1qWJ//2D2mtGz1TpTk8NsumMbNctuDxQkcfX/nVRn67ef+EnvtNq1r51Qt7AbjukpN5++kLmVsXx404qCo96QJ3PrOLZ3f28t6zF/OGlS2jzpHJe/zg8a184NVLp8V0BTOV1fSNOU6pKrmiz2CuSG+mQG+mwH0v7OW7j24lnR96XNZaF6c+4bJ1f/qQNfZKr25r4smt3cPSYq7Dc59966j2cTP9WE3fTAuqStFXPF+JRRzSBQ8/bFIAiEaEhBuZ1BGN/dkC3YMFFjfVHNZ1ip7P+l29rN3azYF0nmQ0wqr5dfi+sn8wT2+mQLbglQNwYzJKLBKhuTZGQ00UVSVf9Cl4StH3yRd94q5DTcwl5jrs68+xuydDwQvyFDyfXNEnF54zVwxes4XgfSbvMZArks57hxXE9/XnOGtJI5ectoCT5teRirmsObHFAvcsNWOD/o+e2M7+wfywNFXFV1AFJXjv+4o/4q8dTxXPUxRwIwIajOjzFXwN8nvhcb4f5A/eB8eM/ONp5K/lqD7GIws/6nil6Cml51KlzyASjHgTAUekXEbXkfLnK31mX4d/Xl8Ptb+0r2J/xT3QintRyuuNs7/yvIejLuGyuClJMhahJhohEXWIuQ6uEwSqygBZ8HwijoR9rId/nvJn8IO0vmyRjXv68BUak1FWzq2l4Omo9ujKr4IdB9L054pAcG+LY3wIEYhFHFJxl95M4Yhq1AA10QjxaPD5ohEhEY0Qd53wJ0Jt3GVOKsiTjEZIxV1S8fA15tJQE6W+xuXAYIFE1OF1J8whFXeJu/ZA0ww3Y4P+tx/ZwqbOgUPmc8KAWfq9UA1G3UUcQYCCr+U8kTBfxBEcERwnSHMkmGvbCd8Do37RRv3aybibo453HRkK9BKULfiCGQqwEUcQETzfDz+TlMvuhMc5YfnL+xyn/PmdsfaLDNsf7BvKO/I6h9rvCEQjDo4j5Is+yVikPMpRRCh6PpmCx/6BPB29GdJ5j3S+yP7BILh7fjAwJxpxgh/XIeoIvipdA3mE0r8PY36GefVxLly9knn1cZ7d0cOOAxkSUaHFjQz7P1Ciqpy9rJHXLJ/D61bMYU4qxmDeY3PnAHHXoSkZoykVHdVbpOD5dKfz9GUKOCJEI6UvLcGNOOSLPpm8R97zaEzGmJOKWXA2x8SMbdPPFoI2zdLvkYRhtTIY2C+ZMWamOm67bIrIPmDbVJdjGmkBRq/SMLvZPRnN7slos+2eLFPVUfPYTPugb4YTkbVjfXvPZnZPRrN7Mprdk4A9vjfGmFnEgr4xxswiFvSPPzdPdQGmIbsno9k9Gc3uCdamb4wxs8q076ff0tKibW1tU10MY4w5Zvb2ZamJRahPHP0cROvWresaq/fOtA/6bW1t2Nw7xpjZYmvXIOd9+SHywH1/fyFNqdhRnUdExuzqbm36xhgzTezty/Kh7wwtKNiYnPjZRqtZLvEk4McVSScAnwEagT8D9oXp16nq3eEx1wJXAR7wKVW9F2OMMew4kOYj332S/QM5/v1PzuZNq+ZOyqwBRx30VXUjcCaAiESAXQSLo/8p8M+q+uXK/CKyGrgMOIVgjdz7RWSVLZlojJnN7vp9B7c9sZ0nth4g6gi3XNnO61eMXm9gokxUm/4FwMuqum2cb6ZLgR+pag7YIiKbgXOB6lc+NsaY49C9z+/hE7c+xfKWFO89ezF//uYTWdhYM6nXnKigfxlwW8X2J0XkQ8Ba4K9VtRtYBDxekWdnmDaKiFwNXA2wdOnSCSqiMcYcG/9094vs6s7wbx88+6B5Ht3cxTU//T2nLqrnpx97fTg1+OSr+kGuiMSAdwI/CZNuAlYQNP10AF8pZR3j8DEHCajqzararqrtra2jehwZY8y0dvPDr3DX+g7W7+zl58/uHravozfDp257mj/+1u9oqInyr5effcwCPkxMTf9i4ClV3QtQegUQkVuAX4SbO4ElFcctBobfDWOMmUHe8fVHAHj7aQso+D4/fHw7X/v1JrIFj0+efyKffPOJJKLHLuDDxAT9y6lo2hGRBaraEW6+G3gufP9z4FYRuZHgQe5K4AmMMWaG+9Yjr/CDx7ex40CGNSfO4XOXnsqK1topKUtVQV9EksCFwEcrkr8kImcSNN1sLe1T1edF5HbgBaAIfMJ67hhjZpLvPzZ8UfqSf7p7A6csrOd7HzmNN65smdIFnKoK+qqaBuaMSLtinPzXA9dXc01jjJmOntx6gM/c+fyY+x695s2T3ivncFVb098K9BMMtiqqaruINBMM2mojqOm/P+y9Y4OzjDEzRvdgnvW7enngxb08vKmLLV2DtNTG6RrIDcv31cvOnDYBHyamTf98Va1cguwa4AFVvUFErgm3P22Ds4wxM4Gqct8Le7n6B+sAiEUc3riqlfe1L+Ydpy9EFb78q43lXjtLm5NTWdxRJmPCtUuB88L33wMeAj6NDc4yxhynCp7P5/7vC/x6Qyf7B3NkCz4An7pgJX/y2qXMrUsMy/+1y8/iQ69bxq83dHL64sapKPJBVRv0FfiViCjwTVW9GZhX6r2jqh0iMjfMa4OzjDHHlYFckbt+v5v/+O1WNuzp56JT5rN0TpKugRwrWmv5+HkrDvpQtr2tmfa25mNc4kOrNuivUdXdYWC/T0Q2jJP3iAZnEa5y097ebqu8GGMmTbbgMZAr0pSMkSl4bOjoY922bn69oZOnt/eQ93xWtKb4wrtO5U9eu2yqi1u1anvv7A5fO0XkDoLmmr2lvvoisgDoDLPb4CxjzJTqTRe4a30Hv3yug99s6hq2L+465Ip+efvk+XV8eE0bF66eR/uypintZjmRqplaOQU4qtofvn8L8DmCQVhXAjeEr3eGh9jgLGPMlHl0cxd//K3fAbC4aag3zdtPX0D7siZ2dGdoTsU4aV4dpy9uYG594mCnOq5VU9OfB9wRfvu5wK2qeo+IPAncLiJXAduB94ENzjLGTBxVZcOeftbv6uX0xQ2cPL+eTN4j5jpsP5CmL1NgaXOSplSM/QM5vvPbLfzbgy+Xj3/4f56P48yMmvuRqmY+/VeAM8ZI308w1fJYx9jgLGPMKNmCx4MbOnlmZw/xiMOuniwdvRnefvpCBnNFFCUWcch7Plu6Bnn05f1s258e95wRRzhjcQPP7uzF85X3nbOYz116KjWxYzvXzXRTTfPOEuD7wHzAB25W1a+KyGexlbOMMSOoKrt7s+zpzeCIEHGEHQcyPLSxk3ue20N/rkg0IhR9RcPuG4++vH/UeRqTUc5Y3MjHz1vBGUsa+er9m/jlc3sAWNCQYEVrLe9/9RJe2N3Hbzbt491nLeKPX7OUs5c2HcuPO22J6tF1jgkf0i5Q1adEpA5YB7wLeD8wcJCVs24jeNi7ELgfOOTgrPb2drWF0Y05fqgqA7kiu3oybOjo5/ndvTy/u48XOvroSRdG5a+Nu1x06nwuPXNhecWorfsH2defY9mcJLVxFxEhV/CIug514fZYPF9xhBnz0LUaIrJOVdtHplfTvNNBMF8+4cPcFzlIv/uQDc4y5jihquQ9n1zRJ53z2NWTYWd3mp3dwWtHb5bOvhyd/bly84tq0Ae76Pn4FXXJmOtw8vw6Lj51PqsX1LOkOYlqMOBpUVMNK1prR00vvKK1dtQslLXxQ4eryCxtpz8SEzIiV0TagLOA3wFrsJWzjJn2Ss0t63f28vzuXl7Y3cfmfQN0D+ZJ5z2K/titAHNSMRY0JpjfkOCMJQ3lmrgACEQdh7qEy8LGGlbNq2NFawo3UvV6TWaCVB30RaQW+Cnwl6raJyI3AZ8n+NL/PMHKWR/BBmcZc0x4vrJ/MEdXf56ugRy9mQK+KgcG8+zszrDjQFBj39Gdpj9bBIIa8omttZy6qIHW2jipeIRkzCXuOiSiERY11rC4qYZFTTUkY5Mxe4s5VqqdZTNKEPB/qKo/A1s5y5gjkS/6bNzTT1+2QMx1mJOK4SvsH8jRny2SiEbwVdnVkyGT96iJRXAdoSddoDudZzBXpC9bZF9/jn39OboGchxI5znYo7qaaIQlzTUsbkpyzrImVs0LAv2rFtQf8xWczNSopveOAN8GXlTVGyvSbeUsMyFKnQxEhGzBI5P3yBV9fFUak1EGckXSueDhXiziEHMdMnmPRNShJhahsy+H4wiuI4hANu8zmC+yvCU1LMAVPJ89vVnqEi57+4Ka8fz6BLmix2DeI50rsqlzgL5MgUQ0QukZYdFXhKD9uSkVQ8Pa9GC+iOs4ZArBsRv39vPk1m5292RIxoIadCru4ghs2jtA3vPH+PSH5jpCKu5Sl3BprYuzbE6Sc9qaaKmN01obo6U2TktdnMaaKI4jNNZEaU7F7CHnLFdNTX8NcAWwXkSeCdOuAy6fDitnbe7sL8+EB8Evtq+K54Oviu8rvoKOaGEKWybLv9h5zydX8MkVg4CTK/pEwp1F3ydb8MrHSdhrQCqOL9W4SlcZ2VtqaP9Q2Yqe4vk+StgmFp5TNShtIhrU9hwRHEfwfL98LFD+bH6YX1WD90rwmcP0ynugpfzha9Ccq/j+0HlK+z1fKXg+BW/EvZPSPecGN+oAACAASURBVAxkiz696Tw1sQhFLzjHgcE8/dkCLbVx0nmP/mwBLyyH5yueKqpaThOBhBshU5i4/yoRR1jcVENNNEJPusD+wdyozzKR6hMu5yxr4o0rW8kWgy+CgZxHwfNZc2ILZyxupKU2Rrbos38gR8QRWmrj1MZd8p6PKkGzSjRCtuhRKCqNqei4vViMOZij7rJ5rBxtl80Lb/x/bOocmIQSHTsiDPszfeQXyZGcxxEpd2UTRmyH+yvzQfDqjNgvAhERohEHN+IEX0aUyjW8YHHXoSEZIxPWfEWgKRmjvsZlX3+OZMyloSZKJPwCizjgOEJEhr7QVJVM3qOhJkoq7pKIRnAE9g/mqUu41MZdCp5PPvxCTkQj5Is+/dki8xviCEG/b0+VRNg+/dLefrZ0DZIt+DSnorTUxlnanGQgV6Q+EaW1Ls6+/hzxqENt3KUmGmF5a4rmVIx8cahnSjQiFDzl5X0D9GYKOCLMScVIxV2KXlCWmliE5mRs1o7+NFNnwrtsTnf/+M5TGMgVy9tuRHAdJwwmlAeHVP4qDgWv0qviRhwSUYe4GyHuOsSjDr4/FEhLo/u0XKsO32tlzXd4FXhkjbhUW4uEZXIdGRUkVLWcL1/0Kfp+uXZcOgaGB3Prrzzx4u7odm8b9GOOJ9O+pi8i+4BtU12OaaQF6DpkrtnJ7s3Y7L6Mbabfl2Wq2joycdoHfTOciKwd6082Y/fmYOy+jG223hcbMWGMMbOIBX1jjJlFLOgff26e6gJMY3Zvxmb3ZWyz8r5Ym74xxswi077LZktLi7a1tU11MYwx5pjJF30cEdzI0Xe5XrduXddYvXemfdBva2vD5tM3xswWv9/Zwzu//lsAtt7wtqM+j4iM2dV92gd9Y4yZLW55+BWuv/vF8rbv64SP5q5mwrWTgB9XJJ0AfAZoxJZLNMaYI/Jn31/LfS+UJylmyxcvmZQR9dWsnLUROBNARCLALuAO4E+Bfz7IcomXAacQLpcoIodcLtEYY2aytmvuGpVWTbPOoUxUl80LgJdVdbzpEsrLJarqFqC0XKIxxsxKIwP+havnTWrAh4lr07+MYNHzElsu0RhjDkJVWX7t3eXtay8+mY++acUxuXbVNX0RiQHvBH4SJt0ErCBo+ukgWC4RjnC5RFVtV9X21tZRPY6MMWZaU1Ue2tg5arpxgKu+++SwgA8cs4APE9O8czHwVGmZRFXdq6qeqvrALQw14dhyicaYWeEna3fy4f94klt+80o5LV/0abvmLh7Y0FlOWzm3dtKbc0aaiOady6lo2rHlEo0xs5mqcv7JcwFwHWfMB7Uweb1zDqXahdGTwIWESyKGvjQdlks0xphjrRTg33b6An5w1bn0ZYqj8nz49W189p2nHOuilVUV9FU1DcwZkXbFOPmvB66v5prGGDMd3fPcnvL7u37fwdcvP6u8feqien7wkdfQlIpNRdGGqbamvxXoJxhsVVTVdhFpJhi01UZQ039/2HvHBmcZY2aMx17ez+W3PH7Q/X/2/bVcfOoCXvmnS6bVGskT8SD3fFU9s2IFmmuAB1R1JfBAuD1ycNZFwDfCQV3GGHPcyBY8tnQNHjTgP/w/z+eNq1rZuLef6+5YP60CPkzO3DuXAueF778HPAR8morBWcAWESkNznpsEspgjDETxvOVFdfdfdD9I3vgfP8j5zKYK3JgMD/ZRTti1QZ9BX4lIgp8U1VvBuaVeu+oaoeIzA3z2uAsY8xxZd22bt5706Nj7mupjfHdPz2XUxc1jLk/FXdJxaffnJbVlmiNqu4OA/t9IrJhnLxHNDiLcFWb9vZ2W+XFGDPhVBURIZ0vsv1AmmXNKV7c08d7vjF2kAe4+1N/wOqF9cewlBOv2t47u8PXThG5g6C5Zm+pr76ILABKIxFscJYxZsp99ufP891Htx7RMZecNp9vfPCcySnQMVbN1MopwFHV/vD9W4DPEQzCuhK4IXy9MzzEBmcZY6bMX/7oaf7rmcOvZ/7jO0/hyte3TV6Bpkg1Nf15wB3hiDIXuFVV7xGRJ4HbReQqYDvwPrDBWcaYiaOq3PH0Lv7q9me56g3L+fu3r8bzFQF+8Pg2nt/dy7UXv4qmVIzHX9nPZTcP72nzP/5wFZ9884lEplnPmmNh2i+M3t7errZcojEzW0dvhtd98dfH5FrHeq6bqSIi6yq60pcddT99EVkiIg+KyIsi8ryI/EWY/lkR2SUiz4Q/l1Qcc62IbBaRjSLy1qO9tjHm+KSq7O3L8sLuPjJ5j7/68TO0XXPXMQn4b1zVOmsC/niqad4pEsyV/5SI1AHrROS+cJ+tnGXMLFf0fP7+zue47YkdAJyysJ5t+9MM5EbPR1Py9tMXsKAhwS2/2TLuuX/x528Y1VVSVTnnC/fz8fNW8N/+4ASKno8bmah1omaOapZL7CCYL5/wYe6LHKTffcgGZxlznNrc2c8f3vhwVedoSsZ4dVvzmD1nNn7hIuLu0AD9v33b6iM+v4jw1N9fWN62gD+2CRk5ICJtwFnA74A12MpZxhwXVJWd3Rk27unnv31/cp6dLWxIcP9fv4lkLAg3UznDpJmAoC8itcBPgb9U1T4RuQn4PMHAq88TrJz1EWxwljHHnKrSnS6wqzvDT5/ayXcf3co5y5o4aX4dG/f089KefvrHaW45mKf//kIak9EpmQ/eVKfaWTajBAH/h6r6MwhWzqrYfwvwi3DTBmcZc5gyeQ8RiLsORV/p6MnSnytQn4gSdx260wUODObp7M/yjQdfZuPefgDOXNLI7p4Mnf05AJKxCOn88Mdm67Z18/K+AVbNreNdZy3i5AV1nDy/juUttTRPg6l/zeSqZnCWAN8GXlTVGyvSbeUsc9zyfUVhWP/toufTny2SLngMZIsoSr7o05MuEI04LG6qYWFjDRFH8H0lXfDIhIG2L1vgiS0HePyV/XSnC0QdIR51iLsRXEd4aW8/L3T04Yd/zzbURMkVPAbD45OxCLmij+cf3h+8z+zoGbZ92auXsriphoaaKBv29LG0OcnFpy2gpTZe/c0yx6VqavprgCuA9SLyTJh2HXD58bJyVmmMwnh/onq+UvB8cgWfdKFIJu8RcaQcFDxfUQVfg2AxNO5BgGCfEu5Xwu0x3ofHBq9Qavka2jd07mTMLV9fBBwRRIIrigier3i+Eo0EeZRgfc5c0afo+Xiq+D7Ba3iB0nkijiCAr0P7NczvjzOmI+/55fuQLXik88HPQK5I0fOD8/la7rmxsLGGzr4sfdkC/dkihfD4aMSh6Puk8x6OCKl4hFTMpTudJ1PwyBV8Cp5PKu4iAt2DBRwHYhEHR2TYfcwVfOKuQ8QRPFX6s8G/34KGBCcvqGdOKsZArkjXQI6tXYP0Z4sM5ot4vjInFccR6M0UyB5G0I1GBNdxyBTG/i89vz7BvPo4RV/L/xb5os/ipho+/Po24m4EJWiKibsOrXVxVGH/QJ5kLMLS5iT1NS59mSI5z6cpGaU5FaM5FcPzleUtqXKbuTHjmbGDs/7opkd5aW8/voYPE4IYXA4KBU/Je/5Q0IQwcEo5gJcCtTl6riNEI0HgFQlqrr7Cvv4cc8KglYq7xNygp0XR84k4QjLm4qsGXxzZIg01UWoTLrGIQ9R1GMgWUKA5GSt/qfmq5X9DwqaRgqcUPR/HEWpjLjWxCFv3DwZt2dkidQmXplSMtpYUTckoCTeCG3E4MJij6CtNyRiJqENLbZxkLEIq7hKR4DM1JqNkCz47utNs25/G832SMZdUPEIiGkFEiLsOZy1p5MS5tdb+bY6pgw3OmrFVgzetauXURQ044S9aZUAQgZjrEI04QS1WdViN2qmoPbuOgxsREtEIyViERNRBFYp+cEDEERxn6LwipZqylt8Pr4kP5ZHKL5uKYFW5DxhWblXIFLyw5jn6LwlfFddxcAQKYe1UCAJg6TOXavQRZ+g6StC04YfniDhScR+ESPgZDiYI7MHVamIRktEINbEg+I3F83VWDoE3ZqpN+5q+iOwDtk11OaaRFqBrqgsxTdm9GZvdl4Obyfdmmaq2jkyc9kHfDCcia8f6k83YvTkYuy8HNxvvjQ1ZM8aYWcSCvjHGzCIW9I8/N091AaYxuzdjs/tycLPu3libvjHGzCLTvstmS0uLtrW1TXUxjDHmuNHS0sK99957r6peNHLftA/6bW1t2MpZxpiZqmsgxwdv+R3/fsU5LG9JTdh5RaRlrPRqVs46qWJ1rGdEpE9E/tJWzjLGmMP3i2d3s3FvP995ZPyFYybKUQd9Vd2oqmeq6pnAOUAauCPc/c+lfap6N4xaOesi4BsiMvZwTWOMmSWK4cj5kZP8rdvWPSnXm6jeOxcAL6vqeCNnyytnqeoWoLRyljHGzFqlyfzccJbWX/x+N2/554e5/ObH2dObnfDrTVTQvwy4rWL7kyLyexH5jog0hWmLgB0VecZdOUtE1orI2n379k1QEY0xZvop1fR39WR4y788zCdvfZpoxOHfPng28+onfgrsqoO+iMSAdwI/CZNuAlYAZxKsofuVUtYxDj/oylmq2q6q7a2to6aOMMaYGaPg+QD88rk9bO4c4KuXncndf/EHXLh63qTMzDoRNf2LgadKK2ap6l5V9VTVB25hqAnHVs4yxpgR/uX+TeX37zlrEZeeuWhSZ6CdiKB/ORVNOyKyoGLfyJWzLhORuIgsx1bOMsbMcg+8WF5dlmVzktzw3tMn/ZrVrpGbBC4kXB0r9KXjZeUsY4yZCqrKbzfv569/8ixz6+Lc8qF2Vi+sJxqZ/Jlxqgr6qpoG5oxIu2Kc/NcD11dzTWOMOV5l8h4/f3YX//HbrWzY08/S5iQ/uOpcls2ZuEFZhzLtR+QaY8zxrCed585ndvPgxk7Wbu1mIFfk5Pl13PCe07j0zEXUxI7tcKVqm3e2Av2ABxRVtV1EmoEfA20EzTvvV9XuMP+1wFVh/k+p6r3VXN8YY6aTHQfS3Pv8Hl7eN8CWrkG2dqXZ0xf0tV/RmuIdZyzkXWcu5NzlzVO2ZvJE1PTPV9XK5cauAR5Q1RtE5Jpw+9MjRuQuBO4XkVXWrm+MOV5t2tvP09t7cBzhwQ2d/PK5DnyFOakYbS0pXn/iHNrmpHjzyXM5dVHDVBcXmJzmnUuB88L33wMeAj5NxYhcYIuIlEbkPjYJZTDGmEmzqyfDV361kTue3kVpdvq6hMufvfEEPvS6NhY11kxtAcdRbdBX4FciosA3VfVmYJ6qdgCoaoeIzA3zLgIerzh23BG5wNUAS5curbKIxhhTPVXlN5u6+K9ndvGL33cAcPUfnMAHXr0ER4QFjQni7vSfTqzaoL9GVXeHgf0+EdkwTt4jGpFLuKJNe3u7rfJijJkyBc/noY37+PqDm3l2Rw/1CZf3nr2IT7555bSu0R9MtV02d4evnSJyB0FzzV4RWRDW8hcAnWF2G5FrjJnWsgWPTXsHeKGjl3XbutncOcCGPf2k8x6Lm2q44T2n8Z6zFxNzj9+VZo866ItICnBUtT98/xbgcwQjb68Ebghf7wwP+Tlwq4jcSPAg10bkGmOmlKpy/4ud/PzZ3by0p5/N+wbKs17WJ1xOWdjAH52zmDec2ML5J889JoOnJls1Nf15wB1htyMXuFVV7xGRJ4HbReQqYDvwPrARucaYief5ykC2SF+2wHV3rOc3m4Y6Ev7d215FQ02UmOvw5V9tZMeBDNdcfDKLGmuYUxvj3x7czG837wegtS7OaYsauHD1PFYvrOdVC+pZ1pzEmcQ5cKbKtF8Yvb29XW25RGNmr3S+yCv7BnluVy8b9vSzcU8/2w+k6csWGMgVqTaE/fFrlvKP7zxlRtTiK4nIOlVtH5leTfPOEuD7wHzAB25W1a+KyGeBPwNKE+FfV7F6lg3OMsYAQdNKX6bIvoEc+/pzdFW87urJ8Mq+QbbuH6Q/Wywfk4xFWDmvjnOXN9NQE6W+Jhq8Jly+99hWntvVx4bPX0TEEQayRfqzRfKeT28mT2Myxge++RjZgs9n3rGad56xEIBEdPr3uJlIR13TDx/SLlDVp0SkDlgHvAt4PzCgql8ekX81wWyc5xIOzgIOOTjLavrGHH9KAb2jL0NHT5ZdPRl2l3+C7X39OfLhXPKVXEeYV5/ghNYUy1tSzG9IsKQpyWmLGlg6Q5tcJsOE1/TDvvil/vj9IvIiB+l3H7LBWcYc53JFj4Fskb19ObYfSLP9wCB7enOk80UG8x77B3Ls6c3S0ZslUxhen3MdYX5DgoWNNby6rYl5DQlaa+O01sVpqXhtrIlaYJ9EEzIiV0TagLOA3wFrCJZL/BCwFvjrcO4dG5xlzDSmqmzuHODBjZ1s2NPPgcE8vZkCvekCvZkC/bki+eLomnkqFiEVd0nFXZqSUV61sJ43nzyX+Q0JFjTUML8hweKmGlpq45O6OIg5PFUHfRGpBX4K/KWq9onITcDnCQZefZ5gucSPYIOzjJmW+rIFntrWzb/+ejPrtnUDML8+QWtdnMZklIWNNWG7eZS6hEtt3GVObYxlzSmWNidpSEan+BOYI1HtLJtRgoD/Q1X9GQTLJVbsvwX4Rbhpg7OMOQaKns9g3qN7MF9+SLp/IEd3usDungw7uzN09GYoeEq+6JdngZxbF+czb1/NxafNZ0HD8TfS1ByeanrvCPBt4EVVvbEifUFp7h1GL5dog7PMjLR/IMfmzgHSeY903mMwXyQbtmkn3AjzGxK4jtCTKSBAKh786u3pzTKQK5KMRRjIFRnMeWH7eJFM3ifigBtxiDoSvEYc0vkifZkC/dki/bkiuaLP3t4sg7ki2aJHwTv4H8cttXEWN9Vw0vw6YuH52lpSrF5Yz2uWN5OM2RIbM101/8JrgCuA9SLyTJh2HXC5LZc4M/m+kvd8ckWffNGn4FW8ej5FTyn6Pvmi4vnBe89XfIWIA5m8z2AuCGhFT4lHHWIRh8G8R9dADt9XohGHuOvgOMJArkh/tkA04lAXd6lNuLiOQ086T182CHa+rzgOgOAIiIAjQiIawQ3bj4NreuSKPkXPpz4RNEcoyr7+HD2ZAi21cZqSUZpTcWIRQUSor4nSUhuj6CkvdvTx25f309mXZW59guZUlMZkjEze4+V9A2zbn56w+xx3HVJxl4Tr4GlwLwueUgjvcTIeoS7hUp+Ikoq71CdcVs5toT4RJRF1SEQjJGMRmpIx5tYHD0fn1MZorIkd19MHmIkxYwdnfeEXL7CnL4sjgq+KKhR9n4KnFH1lrM8tIvjhEGwRUA0CAwTvfS0FM8UPg5mvQ/kFCV+DhCBPmC98D8G3YalMwX7F94MHaRpeqxS8IJjwqVRrjDiCI0FQcmRou3QeLzxPqayVHEfwvGBf1HVwHYeIAxERnIrzqAbnyIeBphRsxupeN1FcR4g4QsHzKRXbdYS6hEvRUwbyQ4NwRKA27hJ3I0Sc0r9NcGd9DUZp5ooexbDGm4q7JGMREtEIEUfoyxTK/75NyRjNqRjd6Tzd6TwHBvPh/4/h5Yu5DucsbWJ5a4o9vVl60nn2D+apiUZYMbeW1QvqOW1RA7UJl1QsuF5NLIIA6bzH7p4MvkJj2P49mCviKyxoSFAbd8kUPFIxl1Q8gjvDBgmZqTHhXTanuy1dg2zZPxjUBCUIxq7jEHWFiOMgUA7QpUCrQCR83KwQ5pHyE2hHhGjEIRENAmTEkfJ5SgF85LkizlCALpWD8LyRcDsSBvEg0IfXDwOZqobXdIIvJR39ReL5lL8ARISIE1yr9KVRoqrBZxfCYB58eXmq5XM5FWWJRhxikeAzu2ENvFQ7j7tB00AsfA1+Sk0QwTGOCK4TnMvzlWQY1FIxFzci5IvBXw3JWISGmigigmrwper5Stx1yqsL+b6SKXgUPJ+6RHTSe4H4vtKbKdA1kMONOCxoSBz1IJ45wJLm5Lh5mo7qzMYcuRkb9L/94VdPdRHMURCR8EtjeLrjSLkd/FhwHKEpFaMpFTtm1zTmWJj2zTsisg/YNtXlmGZagK5D5pq97P4cnN2b8c2U+9MFoKoXjdwx7YO+GU1E1o7VVmcCdn8Ozu7N+GbD/bEnRsYYM4tY0DfGmFnEgv7x6eapLsA0Z/fn4OzejG/G3x9r0zfGmFlk2nfZbGlp0ba2tqkuhjHGTApflT29WerCCe0myrp167pUtXVkejVz75wE/Lgi6QTgM0AjE7hyVltbG7aIijFmpvJ95YTr7ubkxQ3c+ck3TNh5RWTMru7VLKKyETgzPHkE2AXcAfwp8M8HWTnrMuAUwpWzROSQK2cZY8xMVlowZl594thcb4LOcwHwsqqON4iqvHKWqm4BSitnGWPMrPXS3n4AXujoA4KZV9uuuYu2a+4aNX/WRJioBqTLCNa/Lalq5SxjjJktDgzmAfj0RSfTds1dw/Z19mcnfG2Dqmv6IhID3gn8JEy6CVhB0PTTQbByFhzBylkicrWIrBWRtfv27RsrizHGzAjty5pY+3d/yJ/f9vSw9Gc/85ZJWcxmImr6FwNPlVbMmoiVs2y5RGPMbHHDLzfwrUe2lLcf+fT5LG4af1bWakxEm/7lVDTtiMiCin0jV866TETiIrIcWznLGDPLPbn1wLCAD0xqwIfq18hNAhcSro4V+pKtnGWMMePzfeV9//5Yeft/vvUkPnH+iZN+3aqCvqqmCdaIqEy7Ypz81wPXV3NNY4w5nu3qybDmhl8PS9vyxUvKCwZNtmk/ItcYY453m/b2c+V3nuB1K1r46VM7h+178XMXHbOAD9U372wF+glG2BZVtV1EmglG6rYRNO+8P+yyeVQjco0x5ngxmCvyg8e38b1Ht7JyXh0PvzS892Ep4LfUxvjN/3ozNbGjW4KzGhNR0z9fVStXmrkGeEBVbxCRa8LtT9uIXGPMTKOq/NXtz/K6FXP4X//n98P2dfRmh22/aVUr//Se01jYkDimNfuRJqN551LgvPD994CHgE9TMSIX2CIipRG5j41xDmOMmdYqB1Ld8fSuUfv/86rXcO7yZmLu9JrBvtqgr8CvRESBb4b96+epageAqnaIyNww72GPyBWRq4GrAZYuXVplEY0xZmIcGMxz9ufvO+j+b3zwbC45bcFB908H1Qb9Naq6Owzs94nIhnHyHvaIXBucZYyZTroGcrR/4f5R6e85axE3fuDMKSjR0au2y+bu8LVTRO4gaK7ZKyILwlr+AqAzzH7YI3KNMWaqpPNFfrOpi4/+YN1B83zxPadx+bnHZytENfPppwBHVfvD928BPkcw8vZK4Ibw9c7wkJ8Dt4rIjQQPcm1ErjFmyt29voOP//CpQ+Z7xxkL+cr7zph2bfRHqpqa/jzgjvAptAvcqqr3iMiTwO0ichWwHXgf2IhcY8zk8XylO50f1gTz+hVzeOcZCxnIFfnCXS8C8OHXt3HeSa284cQW7nthLx87RLD/h3es5vTFDZy9tGlKe9xMpGm/Rm57e7vaylnGmF+u7zhkkK7G+9sX87dvW00qFsGNHN+1eQARWaeq7SPTq2neWQJ8H5gP+MDNqvpVEfksE7hcojFmdlBVXto7wLcfeYXb1+489AGH4YTWFK/sGzzo/tq4y399Yg0nzq2dkOsdD6pp3ikSLJDylIjUAetEpNSXyZZLNMYAkC14XH7L4zy9vYc3nNjCI5u7Dn3QOP5gZQvfurKduHv0o1lVdcY01xypatbI7SBYJIXwYe6LjL8Slg3OMmaG8H3l+d19vOPrjxzRcYcb8F93why+8v4zaE7FiLvOhAfo2RrwYYJG5IpIG3AW8DtgDbZcojHHJd9XPvPz5/jPx7dPyvlvvuIcTp5fz9I5kztnvDm4qoO+iNQCPwX+UlX7ROQm4PMEA68+T7Bc4kc4wuUSsRG5xhwznq+suO7uCTnXtz7UzptPnovjzN7a9HRW7SybUYKA/0NV/RnYconGTDfZgse2/Wk27OmjPhHl9rU7+OVze7hw9Tx60wW27B9kX39u1HGtdXF+dPVrWdqcxA0D+GxuFpkpqum9I8C3gRdV9caK9AWluXcYvVyiDc4yM5aqUvSVoqcM5IpkCx4RR0jGItQloggwmC9SEw26BBY8n237B3Edh5a6ON2DeTr7s+zpzTGYLxJ3HaIRB1UoeD6NySh1CZeip3SnCxwYzNOXLeD5Sn+2SE86z4+e3AHAyfPr6EkX6E7nyRX9Mcu7bls3K+fWcs7SJu55fg8AD/7NeSxvSR2rW2amQDU1/TXAFcB6EXkmTLsOuNyWSzQQBEHV4D9C0ffp7AtqkyJQGh6SjEUYyBU5MJgnk/dwIw5x18FXZU9vFgWCSqbgSFDTVFV8VXwNzlMKrpWV0FKNVADXERToSRfIFjzm1seZk4rTNZBjV0+GdK5IcypGMu6SjEVIuBGSsQj1NVEaaqKk4i7b9g/yxJYDvNjRR23CpS4RpTkZQ1E2dPTzQkcfmzsHKPpj/2HqhuUreMH+mOugquXtakUjQmMyVt5e2pzk9MVRGpMxGpNR9g/k6ezP8cHXLGXl3FqaUzGrtc9SM3Zw1j3PddCTLpS3PVU8f+in5FAf31dFS686PJANS6vcDt+jQ2m+gjJ2/lJxHAmOLV9Dh64vgBuRcnr5XAxtBx+o4hzlfUPbVBwLY5+LYdvDz+P5PrmCT97zyReD10zeYyBbJFf0UYJg/P+3d7+hkdR3HMffn9nd7N5lk7uo8S6ivSuCniLVKmhpoQ+KCKIgFAon7SPPB4KI4APBR+ozQUEfFRERlKJiS31yFGtRwT4oiH/uFKlcC1oolJ5VzsuZ5JJsvj6YSXZuzSXbbG4zs/N5wZDJzu83+5svs9/5Mcz8fqvHVgb5i9BmLtvTYn6pw+zC8lqC3zfZ5ND+Sa6ZmaTdrFFLEtrNGs1Gjc5KMLfY4aszZwlganeDhaUVvl1cJpG4al+bzko6oNdFu8eYnmyyf7JFu1nn7PLK2rlar4lTc4vMLizTaii6BwAABQtJREFUqCVMtOpMTzSZbDWoJbogT7hYuW37y1lF99SbJ/jnyTND+77VXqiARAJln+V6qCJNMFL3s9VeLHSfHc6qk+TW095ydPdBd9tarzb7za+3ffVr8v+vlUe5uvl2nrsfJOpZgmm36ozV0tsP480a7WaDsXqCBLWs3d32i1oCl060ugHLvm/u7DLtVoOLxhvsHqvTWQnOLneIgP17WtSTJOvVdy+EEuf07Fv1WnqRJX/hSv+LSC/4EbBnV4NdjRr/nV3gf7OLTE80mdnbYlejxun5JeYWO8wvdVhY6jC32OGb+SVOzy8xu7DMzJ4WNx2cOucYTs0tEgFT490etlnRjWzS/92RW1iJWOvF1RORJFr7m+8Tna+HFBEkUjd55ZJg0pPErTymxsfS98hzLm43ufj/3E/+dopZWRT+9o6kL4F/7XQ7CuYSYLDXGkeb47Mxx2djoxKfAxEx3fth4ZO+fZ+k99e7V2cpx2djjs/GRj0+5R9KzszM+uakb2ZWIU765fTcTjeg4ByfjTk+Gxvp+PievplZhbinb2ZWIU76ZmYV4qRfEJKukPSOpL9L+lTSg+uUmZL0uqSPJb0n6bp+65bdgPFpZf8fz+o+PvwjuLAGiU9ue03SR5KO9tYtu0HjI+kLSZ9IOiap3JN2p+O8eNnpBZgBbszWJ4ATwLU9ZZ4EHs3WDwFv9Vu37MuA8RHQztYbpJP9/GSnj6ko8cltfwh4GTi608dTtPiQDh55yU4fx3Ys7ukXRET8JyI+zNZngfWmn7wWeCsr8xlwUNK+PuuW2oDxiYhYHYipkS0j9QTDIPEBkHQ5cAfw/NAaPUSDxmeUOOkXUM/0k3nHgV9mZW4GDpBORtNP3ZGxlfhkty6OASeBv0SE43Pu+fMM8DCw/uD7I2SL8QngTUkfZDP7lZaTfsH0Tj/Zs/kJYCpLXg8AH5HOTdBP3ZGw1fhERCcibiD9Ed/cez97VGwlPpLuBE5GxAfDbe3wDfD7+llE3AjcDtwv6efDavN283P6BZJNP3kU+HPkZiM7T1kBnwM/inRu4r7rltUg8enZ9ijwbUQ8dcEauwO2Gh/gEdIJkZaBFjAJ/DEifnNhWzxc23j+PAacKev5455+QWQn2femn+wps1fS6ni+9wLvZgl/07plN2B8piXtzcrsAm4FPhtGu4dlkPhExCMRcXlEHAQOA2+PYMIf5PwZlzSRlRkHbqM7DWzpjOx4+iV0vuknfwAQEc8C1wAvSeqQTjt5ZKO6EfGnYTV+CAaJzwzwoqQaaUfntYgYtccSB4lPFQwSn33A6+l1gzrwckS8McS2byvf3jEzqxDf3jEzqxAnfTOzCnHSNzOrECd9M7MKcdI3MysQSS9IOilp08dCJT2dDQJ3TNIJSac2reOnd8zMiiN72/cM8FJE9P3muKQHgB9HxD0blXNP38ysQCLiXeDr/GeSrpT0Rjb2z18lHVqn6t3AK5vt3y9nmZkV33PAfRHxD0m3AL8FfrG6UdIB4IfA25vtyEnfzKzAskHifgr8PnsrGKDZU+ww8IeI6Gy2Pyd9M7NiS4BT2Six53MYuL/fnZmZWUFlo3x+LulXkA4eJ+n61e2SrgamgL/1sz8nfTOzApH0CmkCv1rSvyUdAX4NHJF0HPgUuCtX5W7g1ejzUUw/smlmViHu6ZuZVYiTvplZhTjpm5lViJO+mVmFOOmbmVWIk76ZWYU46ZuZVch3x6D3MLcL+HoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Average together all cycles for each beam\n",
    "h_li_avg = {}\n",
    "\n",
    "avgSurfFig, avgSurfAx = plt.subplots(len(beams),1)\n",
    "avgSurfAxCount = 0\n",
    "for beam in beams:\n",
    "    \n",
    "    # if cycles for a given beam are not of the same length, trim ends\n",
    "    cycle_length = np.zeros(len(cycles))\n",
    "    for cycleInd in np.arange(0, len(cycles)):\n",
    "        cycle_length[cycleInd] = len(x_atc[cycles[cycleInd]][beam])\n",
    "    \n",
    "    if np.any(cycle_length != cycle_length[0]):\n",
    "        print('oops! Not all cycles in beam {} are the same length. Trimming'.format(beam))\n",
    "        long_cycle = int(np.where(cycle_length == np.max(cycle_length))[0]) # identify which is longer. np.where returns a tuple of arrays, so index 0\n",
    "        short_cycle = int(np.where(cycle_length != np.max(cycle_length))[0])\n",
    "        if x_atc[cycles[long_cycle]][beam][0] is not x_atc[cycles[short_cycle]][beam][0]: #trim beginning of array\n",
    "            trimmed_x_atc = x_atc[cycles[long_cycle]][beam][1:]\n",
    "            trimme_h_li = h_li[cycles[long_cycle]][beam][1:]\n",
    "        \n",
    "        if x_atc[cycles[long_cycle]][beam][-1] is not x_atc[cycles[short_cycle]][beam][-1]: #trim end of array\n",
    "            trimmed_x_atc = x_atc[cycles[long_cycle]][beam][:-1]\n",
    "            trimmed_h_li = h_li[cycles[long_cycle]][beam][:-1]\n",
    "            \n",
    "        h_li[cycles[long_cycle]][beam] = trimmed_h_li\n",
    "        x_atc[cycles[long_cycle]][beam] = trimmed_x_atc\n",
    "        \n",
    "        \n",
    "    # average h_li across all cycles for the beam\n",
    "    cycleSum = np.zeros(np.shape(x_atc[cycles[0]][beam]))\n",
    "    cycleCount = 0.\n",
    "    for cycle in cycles:\n",
    "        cycleSum += h_li[cycle][beam]\n",
    "        cycleCount += 1.\n",
    "    \n",
    "    h_li_avg[beam] = cycleSum / cycleCount \n",
    "    \n",
    "    # plot up average surface for each beam\n",
    "    avgSurfAx[avgSurfAxCount].plot(x_atc[cycles[0]][beam], h_li_avg[beam]) \n",
    "    avgSurfAxCount += 1\n",
    "    \n",
    "    #subtract average surface off all cycles for each beam\n",
    "    for cycle in cycles:\n",
    "        h_li[cycle][beam] -= h_li_avg[beam]\n",
    "          "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "x and y must have same first dimension, but have shapes (21161,) and (21615,)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-114-1f058366859b>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx_atc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mcycles\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mbeam\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mh_li\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mcycles\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mbeams\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshow\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/srv/conda/envs/notebook/lib/python3.7/site-packages/matplotlib/pyplot.py\u001b[0m in \u001b[0;36mplot\u001b[0;34m(scalex, scaley, data, *args, **kwargs)\u001b[0m\n\u001b[1;32m   2761\u001b[0m     return gca().plot(\n\u001b[1;32m   2762\u001b[0m         *args, scalex=scalex, scaley=scaley, **({\"data\": data} if data\n\u001b[0;32m-> 2763\u001b[0;31m         is not None else {}), **kwargs)\n\u001b[0m\u001b[1;32m   2764\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2765\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/srv/conda/envs/notebook/lib/python3.7/site-packages/matplotlib/axes/_axes.py\u001b[0m in \u001b[0;36mplot\u001b[0;34m(self, scalex, scaley, data, *args, **kwargs)\u001b[0m\n\u001b[1;32m   1644\u001b[0m         \"\"\"\n\u001b[1;32m   1645\u001b[0m         \u001b[0mkwargs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcbook\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnormalize_kwargs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmlines\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mLine2D\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1646\u001b[0;31m         \u001b[0mlines\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_get_lines\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1647\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0mline\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mlines\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1648\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd_line\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mline\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/srv/conda/envs/notebook/lib/python3.7/site-packages/matplotlib/axes/_base.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m    214\u001b[0m                 \u001b[0mthis\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    215\u001b[0m                 \u001b[0margs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 216\u001b[0;31m             \u001b[0;32myield\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_plot_args\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mthis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    217\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    218\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mget_next_color\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/srv/conda/envs/notebook/lib/python3.7/site-packages/matplotlib/axes/_base.py\u001b[0m in \u001b[0;36m_plot_args\u001b[0;34m(self, tup, kwargs)\u001b[0m\n\u001b[1;32m    340\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    341\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 342\u001b[0;31m             raise ValueError(f\"x and y must have same first dimension, but \"\n\u001b[0m\u001b[1;32m    343\u001b[0m                              f\"have shapes {x.shape} and {y.shape}\")\n\u001b[1;32m    344\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndim\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m2\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndim\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: x and y must have same first dimension, but have shapes (21161,) and (21615,)"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAANQklEQVR4nO3cX2id933H8fdndg3rnzWhUUtnp9QbTlNfNCNR0zDWLV3ZamcXptCLpKVhoWDCmtLLhMHai9ysF4NSktSYYEJv6os1tO5IGwajzSBLFxlSJ05I0VwWay7EaUsHKSw4+e7inE1Cka3H5xxJjr7vFwj0nOcn6asf8tuPj3WeVBWSpO3vd7Z6AEnS5jD4ktSEwZekJgy+JDVh8CWpCYMvSU2sG/wkx5K8nOS5i5xPkm8kWUxyKsmNsx9TkjStIVf4jwAHLnH+ILBv/HYY+Ob0Y0mSZm3d4FfVE8CvLrHkEPCtGnkKuCrJ+2c1oCRpNnbO4HPsBs6uOF4aP/aL1QuTHGb0rwDe8Y533HT99dfP4MtLUh8nT558parmJvnYWQQ/azy25v0aquoocBRgfn6+FhYWZvDlJamPJP856cfO4rd0loBrVxzvAc7N4PNKkmZoFsE/Adw5/m2dW4DfVNWbns6RJG2tdZ/SSfJt4FbgmiRLwFeBtwFU1RHgMeA2YBH4LXDXRg0rSZrcusGvqjvWOV/AF2c2kSRpQ/hKW0lqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpoYFPwkB5K8mGQxyX1rnH93ku8n+WmS00numv2okqRprBv8JDuAB4GDwH7gjiT7Vy37IvB8Vd0A3Ar8Q5JdM55VkjSFIVf4NwOLVXWmql4DjgOHVq0p4F1JArwT+BVwYaaTSpKmMiT4u4GzK46Xxo+t9ADwYeAc8Czw5ap6Y/UnSnI4yUKShfPnz084siRpEkOCnzUeq1XHnwKeAX4f+CPggSS/96YPqjpaVfNVNT83N3fZw0qSJjck+EvAtSuO9zC6kl/pLuDRGlkEfg5cP5sRJUmzMCT4TwP7kuwd/0fs7cCJVWteAj4JkOR9wIeAM7McVJI0nZ3rLaiqC0nuAR4HdgDHqup0krvH548A9wOPJHmW0VNA91bVKxs4tyTpMq0bfICqegx4bNVjR1a8fw74y9mOJkmaJV9pK0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqYlDwkxxI8mKSxST3XWTNrUmeSXI6yY9nO6YkaVo711uQZAfwIPAXwBLwdJITVfX8ijVXAQ8BB6rqpSTv3aiBJUmTGXKFfzOwWFVnquo14DhwaNWazwKPVtVLAFX18mzHlCRNa0jwdwNnVxwvjR9b6Trg6iQ/SnIyyZ1rfaIkh5MsJFk4f/78ZBNLkiYyJPhZ47FadbwTuAn4K+BTwN8lue5NH1R1tKrmq2p+bm7usoeVJE1u3efwGV3RX7vieA9wbo01r1TVq8CrSZ4AbgB+NpMpJUlTG3KF/zSwL8neJLuA24ETq9Z8D/h4kp1J3g58DHhhtqNKkqax7hV+VV1Icg/wOLADOFZVp5PcPT5/pKpeSPJD4BTwBvBwVT23kYNLki5PqlY/Hb855ufna2FhYUu+tiS9VSU5WVXzk3ysr7SVpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpiUHBT3IgyYtJFpPcd4l1H03yepLPzG5ESdIsrBv8JDuAB4GDwH7gjiT7L7Lua8Djsx5SkjS9IVf4NwOLVXWmql4DjgOH1lj3JeA7wMsznE+SNCNDgr8bOLvieGn82P9Lshv4NHDkUp8oyeEkC0kWzp8/f7mzSpKmMCT4WeOxWnX8deDeqnr9Up+oqo5W1XxVzc/NzQ2dUZI0AzsHrFkCrl1xvAc4t2rNPHA8CcA1wG1JLlTVd2cypSRpakOC/zSwL8le4L+A24HPrlxQVXv/7/0kjwD/ZOwl6cqybvCr6kKSexj99s0O4FhVnU5y9/j8JZ+3lyRdGYZc4VNVjwGPrXpszdBX1V9PP5YkadZ8pa0kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqYlBwU9yIMmLSRaT3LfG+c8lOTV+ezLJDbMfVZI0jXWDn2QH8CBwENgP3JFk/6plPwf+rKo+AtwPHJ31oJKk6Qy5wr8ZWKyqM1X1GnAcOLRyQVU9WVW/Hh8+BeyZ7ZiSpGkNCf5u4OyK46XxYxfzBeAHa51IcjjJQpKF8+fPD59SkjS1IcHPGo/VmguTTzAK/r1rna+qo1U1X1Xzc3Nzw6eUJE1t54A1S8C1K473AOdWL0ryEeBh4GBV/XI240mSZmXIFf7TwL4ke5PsAm4HTqxckOQDwKPA56vqZ7MfU5I0rXWv8KvqQpJ7gMeBHcCxqjqd5O7x+SPAV4D3AA8lAbhQVfMbN7Yk6XKlas2n4zfc/Px8LSwsbMnXlqS3qiQnJ72g9pW2ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNTEo+EkOJHkxyWKS+9Y4nyTfGJ8/leTG2Y8qSZrGusFPsgN4EDgI7AfuSLJ/1bKDwL7x22HgmzOeU5I0pSFX+DcDi1V1pqpeA44Dh1atOQR8q0aeAq5K8v4ZzypJmsLOAWt2A2dXHC8BHxuwZjfwi5WLkhxm9C8AgP9J8txlTbt9XQO8stVDXCHci2XuxTL3YtmHJv3AIcHPGo/VBGuoqqPAUYAkC1U1P+Drb3vuxTL3Ypl7scy9WJZkYdKPHfKUzhJw7YrjPcC5CdZIkrbQkOA/DexLsjfJLuB24MSqNSeAO8e/rXML8Juq+sXqTyRJ2jrrPqVTVReS3AM8DuwAjlXV6SR3j88fAR4DbgMWgd8Cdw342kcnnnr7cS+WuRfL3Itl7sWyifciVW96ql2StA35SltJasLgS1ITGx58b8uwbMBefG68B6eSPJnkhq2YczOstxcr1n00yetJPrOZ822mIXuR5NYkzyQ5neTHmz3jZhnwZ+TdSb6f5KfjvRjy/4VvOUmOJXn5Yq9VmribVbVhb4z+k/c/gD8AdgE/BfavWnMb8ANGv8t/C/CTjZxpq94G7sUfA1eP3z/YeS9WrPsXRr8U8JmtnnsLfy6uAp4HPjA+fu9Wz72Fe/G3wNfG788BvwJ2bfXsG7AXfwrcCDx3kfMTdXOjr/C9LcOydfeiqp6sql+PD59i9HqG7WjIzwXAl4DvAC9v5nCbbMhefBZ4tKpeAqiq7bofQ/aigHclCfBORsG/sLljbryqeoLR93YxE3Vzo4N/sVsuXO6a7eByv88vMPobfDtady+S7AY+DRzZxLm2wpCfi+uAq5P8KMnJJHdu2nSba8hePAB8mNELO58FvlxVb2zOeFeUibo55NYK05jZbRm2gcHfZ5JPMAr+n2zoRFtnyF58Hbi3ql4fXcxtW0P2YidwE/BJ4HeBf0vyVFX9bKOH22RD9uJTwDPAnwN/CPxzkn+tqv/e6OGuMBN1c6OD720Zlg36PpN8BHgYOFhVv9yk2TbbkL2YB46PY38NcFuSC1X13c0ZcdMM/TPySlW9Crya5AngBmC7BX/IXtwF/H2NnsheTPJz4Hrg3zdnxCvGRN3c6Kd0vC3DsnX3IskHgEeBz2/Dq7eV1t2LqtpbVR+sqg8C/wj8zTaMPQz7M/I94ONJdiZ5O6O71b6wyXNuhiF78RKjf+mQ5H2M7hx5ZlOnvDJM1M0NvcKvjbstw1vOwL34CvAe4KHxle2F2oZ3CBy4Fy0M2YuqeiHJD4FTwBvAw1W17W4tPvDn4n7gkSTPMnpa496q2na3TU7ybeBW4JokS8BXgbfBdN301gqS1ISvtJWkJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5Ka+F/Xe3Wlc9XddQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x_atc[cycles[0]][beam], h_li[cycles[0]][beams[3]])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/ipykernel_launcher.py:22: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8ffc8c49e64d45a9b8c159342c1ead72",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/ipykernel_launcher.py:22: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cc045cf1c414439b973f35c1fa7d4234",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/ipykernel_launcher.py:22: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8bc0fe2d6f3c4680b7636c6f80e78c1c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/ipykernel_launcher.py:22: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "407c1ae9446e4a49b694e8cd951e540b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/ipykernel_launcher.py:22: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "13ff755f8bce489ea9af4e67401ff20c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/ipykernel_launcher.py:22: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "35b0556518ad443ba4fca0d3b0773cb4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_veloc = len(cycles) - 1\n",
    "\n",
    "segment_length = 10000 # m\n",
    "x1 = 2.915e7#x_atc[cycles[0]][beams[0]][1000] <-- the very first x value in a file; doesn't work, I think b/c nans # 2.93e7\n",
    "\n",
    "search_width = 2000 # m\n",
    "dx = 20 # meters between x_atc points\n",
    "\n",
    "for veloc_number in range(n_veloc):\n",
    "    cycle1 = cycles[veloc_number]\n",
    "    cycle2 = cycles[veloc_number+1]\n",
    "    t1_string = times[cycle1]['gt1l'][0].astype(str) #figure out later if just picking hte first one it ok\n",
    "    t1 = Time(t1_string)\n",
    "                       \n",
    "    t2_string = times[cycle2]['gt1l'][0].astype(str) #figure out later if just picking hte first one it ok\n",
    "    t2 = Time(t2_string)\n",
    "    \n",
    "    dt = (t2 - t1).jd # difference in julian days\n",
    "        \n",
    "    velocities = {}     \n",
    "    for beam in beams:\n",
    "        fig1, axs = plt.subplots(4,1)\n",
    "        \n",
    "        # cut out small chunk of data at time t1 (first cycle)\n",
    "        x_full_t1 = x_atc[cycle1][beam]\n",
    "        ix_x1 = np.arange(len(x_full_t1))[x_full_t1 >= x1][0]\n",
    "        ix_x2 = ix_x1 + int(np.round(segment_length/dx))      \n",
    "        x_t1 = x_full_t1[ix_x1:ix_x2]\n",
    "        h_li1 = h_li_diff[cycle1][beam][ix_x1-1:ix_x2-1] # start 1 index earlier because the data are differentiated\n",
    "        \n",
    "        # cut out a wider chunk of data at time t2 (second cycle)\n",
    "        x_full_t2 = x_atc[cycle2][beam]\n",
    "        ix_x3 = ix_x1 - int(np.round(search_width/dx)) # offset on earlier end by # indices in search_width\n",
    "        ix_x4 = ix_x2 + int(np.round(search_width/dx)) # offset on later end by # indices in search_width\n",
    "        x_t2 = x_full_t2[ix_x3:ix_x4]\n",
    "        h_li2 = h_li_diff[cycle2][beam][ix_x3:ix_x4]\n",
    "\n",
    "        # plot data\n",
    "        axs[0].plot(x_t2, h_li2, 'r')\n",
    "        axs[0].plot(x_t1, h_li1, 'k')\n",
    "        axs[0].set_xlabel('x_atc (m)')\n",
    "        \n",
    "        # correlate old with newer data\n",
    "        corr = correlate(h_li1, h_li2, mode = 'valid', method = 'direct') \n",
    "        norm_val = np.sqrt(np.sum(h_li1**2)*np.sum(h_li2**2)) # normalize so values range between 0 and 1\n",
    "        corr = corr / norm_val\n",
    "        \n",
    "        \n",
    "#         lagvec = np.arange( -(len(h_li1) - 1), len(h_li2), 1)# for mode = 'full'\n",
    "#         lagvec = np.arange( -int(search_width/dx) - 1, int(search_width/dx) +1, 1) # for mode = 'valid'\n",
    "        lagvec = np.arange(- int(np.round(search_width/dx)), int(search_width/dx) +1,1)# for mode = 'valid'\n",
    "\n",
    "        shift_vec = lagvec * dx\n",
    "\n",
    "        ix_peak = np.arange(len(corr))[corr == np.nanmax(corr)][0]\n",
    "        best_lag = lagvec[ix_peak]\n",
    "        best_shift = shift_vec[ix_peak]\n",
    "        velocities[beam] = best_shift/(dt/365)\n",
    "\n",
    "        axs[1].plot(lagvec,corr)\n",
    "        axs[1].plot(lagvec[ix_peak],corr[ix_peak], 'r*')\n",
    "        axs[1].set_xlabel('lag (samples)')\n",
    "\n",
    "        axs[2].plot(shift_vec,corr)\n",
    "        axs[2].plot(shift_vec[ix_peak],corr[ix_peak], 'r*')\n",
    "        axs[2].set_xlabel('shift (m)')\n",
    "\n",
    "        # plot shifted data\n",
    "        axs[3].plot(x_t2, h_li2, 'r')\n",
    "        axs[3].plot(x_t1 - best_shift, h_li1, 'k')\n",
    "        axs[3].set_xlabel('x_atc (m)')\n",
    "        \n",
    "        axs[0].text(x_t2[100], 0.6*np.nanmax(h_li2), beam)\n",
    "        axs[1].text(lagvec[5], 0.6*np.nanmax(corr), 'best lag: ' + str(best_lag) + '; corr val: ' + str(np.round(corr[ix_peak],3)))\n",
    "        axs[2].text(shift_vec[5], 0.6*np.nanmax(corr), 'best shift: ' + str(best_shift) + ' m'+ '; corr val: ' + str(np.round(corr[ix_peak],3)))\n",
    "        axs[2].text(shift_vec[5], 0.3*np.nanmax(corr), 'veloc of ' + str(np.round(best_shift/(dt/365),1)) + ' m/yr')\n",
    "\n",
    "        \n",
    "    fig1.suptitle('black = older cycle data, red = newer cycle data to search across')\n",
    "\n",
    "\n",
    "        \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:notebook] *",
   "language": "python",
   "name": "conda-env-notebook-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
